#!/usr/bin/python3 """ Joint radio and gamma-ray (jrag) timing analysis. Basically a wrapper around the PINT software. Author: Lars Nieder """ ################################################################################ # IMPORT: Modules ################################################################################ from pint.models import get_model from pint.fermi_toas import load_Fermi_TOAs from pint.observatory.satellite_obs import SatelliteObs from pint.templates.lctemplate import LCTemplate,prim_io from pint.templates.lcprimitives import LCGaussian from pint.residuals import Residuals from pint import toa from astropy import units as u from astropy.coordinates import Angle from astropy.time import Time import emcee from matplotlib.gridspec import GridSpec import matplotlib.pyplot as plt import matplotlib as mpl import corner from optparse import OptionParser,OptionGroup import sys, os, shutil import time import numpy as np ################################################################################ # CONSTANTS: Useful constants ################################################################################ two_sqrt_two_log_two = 2 * np.sqrt(2 * np.log(2)) lob_fwhm = 0.005*two_sqrt_two_log_two # sigma = 0.5% of a rotation hib_fwhm = 0.50*two_sqrt_two_log_two # sigma = 50% of a rotation ################################################################################ # FUNCTIONS: Useful helper functions ################################################################################ def order_of_mag(values): mask = np.where(np.array(values) != 0.0)[0] print(values) print(np.shape(values)) if len(mask): print(mask) print(np.shape(mask)) mags = np.zeros_like(values) # taking care of special case: value = 0 mags[mask] = np.floor(np.log10(np.abs(np.array(values)[mask]))) else: mags = np.floor(np.log10(np.abs(np.array(values)))) return mags.astype(int) ################################################################################ # FUNCTIONS: Analysis plots ################################################################################ def plot_phases(fermi_toas,weights,radio_toas,timing_model,isolated=False,template=None,savepath=None): # Folding radio_residuals = Residuals(radio_toas,timing_model).phase_resids radio_toa_uncerts = (radio_toas.get_errors()*timing_model.F0.quantity).to_value(u.dimensionless_unscaled) ntoa = len(radio_residuals) fermi_phases = np.mod(timing_model.phase(fermi_toas).frac,1.0).value npho = len(fermi_phases) # Computing statistics and relevant quantities rms = np.sqrt(np.sum(np.array(radio_residuals)**2)/ntoa)/timing_model.F0.value chi2 = chi_square(radio_residuals,radio_toa_uncerts) redchi2 = red_chi_square(radio_residuals,radio_toa_uncerts) sumw2 = np.sum(weights**2) ht = htest(fermi_phases,weights) ll = loglike(fermi_phases,weights,template) # Radio: Getting other properties (orbital phase, day-of-year, etc.) radio_t = timing_model.get_barycentric_toas(radio_toas).value rstart = np.min(radio_t) rfinish = np.max(radio_t) t = Time(timing_model.get_barycentric_toas(radio_toas),format='mjd') radio_doy = (t - Time(np.floor(t.decimalyear),format='decimalyear')).to_value(u.d) radio_orbphi = timing_model.orbital_phase(radio_toas,anom="mean",radians=False) # Gamma: Getting other properties (orbital phase, day-of-year, etc.) fermi_orbphi = timing_model.orbital_phase(fermi_toas,anom="mean",radians=False) fermi_t = timing_model.get_barycentric_toas(fermi_toas).value gstart = np.min(fermi_t) gfinish = np.max(fermi_t) t = Time(radio_toas.get_mjds(),format='mjd') fermi_doy = (t - Time(np.floor(t.decimalyear),format='decimalyear')).to_value(u.d) # Plotting fig = plt.figure(figsize=(16,10)) gs = GridSpec(80, 50, figure=fig) if isolated == False: ax0 = fig.add_subplot(gs[19:,0:11]) ax01 = fig.add_subplot(gs[0:18,0:11]) ax1 = fig.add_subplot(gs[19:,14:25]) ax2 = fig.add_subplot(gs[19:36,29:]) ax3 = fig.add_subplot(gs[41:58,29:]) ax4 = fig.add_subplot(gs[63:80,29:]) axtxt = fig.add_subplot(gs[0:18,14:]) else: ax0 = fig.add_subplot(gs[19:,0:11]) ax01 = fig.add_subplot(gs[0:18,0:11]) ax2 = fig.add_subplot(gs[19:47,15:]) ax4 = fig.add_subplot(gs[52:80,15:]) axtxt = fig.add_subplot(gs[0:18,15:]) phi = np.append(fermi_phases, fermi_phases + 1.0) t = np.append(fermi_t,fermi_t) w = np.append(weights,weights) o = np.append(fermi_orbphi,fermi_orbphi) i = np.argsort(w) # Phase histogram xbins=50 wcounts,bins = np.histogram(fermi_phases, xbins, weights=weights) sqerrors,bins = np.histogram(fermi_phases, xbins, weights=(weights ** 2.0)) width = (bins[1] - bins[0]) bg = (weights.sum() - ((weights ** 2.0).sum()))/xbins errors = np.sqrt(1.0 + sqerrors) ax01.bar(np.append(bins[:-1],bins[:-1] + 1.0),np.append(wcounts,wcounts),width=width,color="lightgray",edgecolor="lightgray",capsize=0,linewidth=1.0,align='edge') ax01.errorbar(np.append(bins[:-1],bins[:-1] + 1.0)+width/2.0,np.append(wcounts,wcounts),yerr=np.append(errors,errors),ecolor="k",elinewidth=1.0,color='none') ax01.step(np.append(bins[:-1],bins+ 1.0),np.append(np.append(wcounts,wcounts),wcounts[0]),color='k',where='post',linewidth=1.0) ax01.hlines(bg,0.0,2.0,linestyle="--",color="steelblue",linewidth=1.5) if template!=None: ticks = np.arange(0.0,1.0,0.001) template_curve = bg + (weights ** 2.0).sum()/xbins * template(ticks) ax01.plot(ticks,template_curve,color='coral',linewidth=1.5) ax01.set_xlim(xmin=0.0,xmax=2.0) ax01.set_ylim(ymin=0.0,ymax=(wcounts + errors).max() * 1.1) ax01.set_ylabel('Weighted Counts') ax01.tick_params(axis='both', which='major') ax01.set_xticklabels([]) # Phase and orbital phase vs time plots ax0.scatter(phi[i],t[i],5,c=w[i],vmin=0,vmax=1,cmap='binary') if isolated == False: ax1.scatter(phi[i],o[i],5,c=w[i],vmin=0,vmax=1,cmap='binary') ax0.errorbar(radio_residuals + 1,radio_t,xerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') ax0.set_xlabel('Spin phase') ax0.set_ylabel('Time (MJD)') ax0.set_xlim(0,2) ax0.set_ylim(bottom=np.minimum(fermi_t.min(),radio_t.min()), top=np.maximum(fermi_t.max(),radio_t.max())) if isolated == False: ax1.errorbar(radio_residuals + 1,radio_orbphi,xerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') ax1.axhline(0.25,linestyle="--",color="steelblue",linewidth=1.5) ax1.set_xlabel('Spin phase') ax1.set_ylabel('Orbital phase') ax1.set_xlim(0,2) ax1.set_ylim(0,1) ax2.errorbar(radio_t,radio_residuals,yerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') ax2.axhline(0.0,linestyle="--",color="steelblue",linewidth=1.5) ax2.set_xlabel('Time (MJD)') ax2.set_ylabel('Residuals') if isolated == False: ax3.errorbar(radio_orbphi,radio_residuals,yerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') ax3.axvline(0.25,linestyle="--",color="steelblue",linewidth=1.5) ax3.axhline(0.0,linestyle="--",color="steelblue",linewidth=1.5) ax3.set_xlabel('Orbital phase') ax3.set_ylabel('Residuals') ax4.errorbar(radio_doy,radio_residuals,yerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') ax4.axhline(0.0,linestyle="--",color="steelblue",linewidth=1.5) ax4.set_xlabel('Day of year') ax4.set_ylabel('Residuals') ax4.set_xlim(0,365.25) axtxt.set_xticks([]) axtxt.set_yticks([]) # Some basic pulsar parameters axtxt.text(0.01,0.90,"PSR {0:s}".format(timing_model.PSR.quantity)) axtxt.text(0.01,0.80,"F0 = {0:.2f} Hz".format(timing_model.F0.value)) axtxt.text(0.01,0.70,"F1 = {0:.2e} Hz/s".format(timing_model.F1.value)) if isolated == False: axtxt.text(0.01,0.60,"PB = {0:.2f} d".format(timing_model.PB.value)) # Radio statistics axtxt.text(0.51,0.90,"Radio ({0:d} TOAs):".format(ntoa)) axtxt.text(0.51,0.80,"rms = ${0:.1f} \mu s$".format(rms*1e06)) axtxt.text(0.51,0.70,"$\chi^2 = {0:.1f}$".format(chi2)) axtxt.text(0.51,0.60,"$\chi^2/({0:d}-0) = {1:.1f}$".format(ntoa,redchi2)) axtxt.text(0.51,0.50,"Start MJD {0:.0f}".format(np.floor(rstart))) axtxt.text(0.51,0.40,"Finish MJD {0:.0f}".format(np.ceil(rfinish))) # Gamma statistics axtxt.text(0.76,0.90,"Gamma ({0:d} Photons):".format(npho)) axtxt.text(0.76,0.80,"$\Sigma w^2 = {0:.1f}$".format(sumw2)) axtxt.text(0.76,0.70,"$H = {0:.1f}$".format(ht)) axtxt.text(0.76,0.60,"log $\mathcal{L}"+" = {0:.1f}$".format(ll)) axtxt.text(0.76,0.50,"Start MJD {0:.0f}".format(np.floor(gstart))) axtxt.text(0.76,0.40,"Finish MJD {0:.0f}".format(np.ceil(gfinish))) #gs.tight_layout(fig) if savepath is not None: plt.savefig(savepath,bbox_inches='tight') plt.show() plt.close() def plot_radio_residuals(radio_toas,timing_model,orbphi_range=0.0,error_range=0.0,isolated=False,savepath=None,debug=0): # Folding radio_residuals = Residuals(radio_toas,timing_model).phase_resids radio_toa_uncerts = (radio_toas.get_errors()*timing_model.F0.quantity).to_value(u.dimensionless_unscaled) ntoa = len(radio_residuals) # Computing statistics and relevant quantities chi2 = chi_square(radio_residuals,radio_toa_uncerts) redchi2 = red_chi_square(radio_residuals,radio_toa_uncerts) rms = np.sqrt(np.sum(np.array(radio_residuals)**2)/ntoa)/timing_model.F0.value # Getting other properties (orbital phase, day-of-year, etc.) radio_t = timing_model.get_barycentric_toas(radio_toas).value t = Time(timing_model.get_barycentric_toas(radio_toas),format='mjd') start = np.min(radio_t) tspan = np.max(radio_t) - start radio_doy = (t - Time(np.floor(t.decimalyear),format='decimalyear')).to_value(u.d) radio_orbphi = timing_model.orbital_phase(radio_toas,anom="mean",radians=False) notok = None if error_range: error_notok = radio_toas.table["error"] > error_range if (isolated == False) and orbphi_range: orbphi_notok = (radio_orbphi > 0.25-0.5*orbphi_range) & (radio_orbphi < 0.25+0.5*orbphi_range) notok = orbphi_notok+error_notok else: notok = error_notok else: if (isolated == False) and orbphi_range: orbphi_notok = (radio_orbphi > 0.25-0.5*orbphi_range) & (radio_orbphi < 0.25+0.5*orbphi_range) notok = orbphi_notok # Plotting fig = plt.figure(figsize=(10,10)) gs = GridSpec(56, 56, figure=fig) if isolated == False: ax1 = fig.add_subplot(gs[0:16,:]) ax2 = fig.add_subplot(gs[20:36,:]) ax3 = fig.add_subplot(gs[40:56,:]) else: ax1 = fig.add_subplot(gs[0:26,:]) ax3 = fig.add_subplot(gs[30:56,:]) # Residual phase vs time plot ax1.errorbar(radio_t,radio_residuals,yerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') if notok is not None: ax1.errorbar(radio_t[notok],radio_residuals[notok],\ yerr=radio_toa_uncerts[notok],\ ls='none',color='steelblue',marker='x') ax1.axhline(0.0,linestyle="--",color="steelblue",linewidth=1.5) ax1.set_xlabel('Time (MJD)') ax1.set_ylabel('Residuals') textstr = '\n'.join(( "#TOA = ${0:d}$".format(ntoa), "$\chi^2 = {0:.1f}$".format(chi2), "$\chi^2/({0:d}-0) = {1:.1f}$".format(ntoa,redchi2), "rms = ${0:.1f} \mu s$".format(rms*1e06) )) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='white', alpha=0.8) # place a text box in upper left in axes coords ax1.text(0.02, 0.95, textstr, transform=ax1.transAxes, verticalalignment='top', bbox=props) if isolated == False: # Residual phase vs orbital phase plot ax2.errorbar(radio_orbphi,radio_residuals,yerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') if notok is not None: ax2.errorbar(radio_orbphi[notok],radio_residuals[notok],\ yerr=radio_toa_uncerts[notok],\ ls='none',color='steelblue',marker='x') if orbphi_range: ax2.axvspan(0.25-0.5*orbphi_range,0.25+0.5*orbphi_range,color="steelblue",alpha=0.2) ax2.axvline(0.25,linestyle="--",color="steelblue",linewidth=1.5) ax2.axhline(0.0,linestyle="--",color="steelblue",linewidth=1.5) ax2.set_xlabel('Orbital phase') ax2.set_ylabel('Residuals') # Residual phase vs day of year plot ax3.errorbar(radio_doy,radio_residuals,yerr=radio_toa_uncerts,\ ls='none',color='red',marker='x') if notok is not None: ax3.errorbar(radio_doy[notok],radio_residuals[notok],\ yerr=radio_toa_uncerts[notok],\ ls='none',color='steelblue',marker='x') ax3.axhline(0.0,linestyle="--",color="steelblue",linewidth=1.5) ax3.set_xlabel('Day of year') ax3.set_ylabel('Residuals') ax3.set_xlim(0,365.25) #gs.tight_layout(fig) if savepath is not None: plt.savefig(savepath,bbox_inches='tight') if debug > 1: plt.show() else: plt.show() plt.close() def plot_chains(samples,errors,labels,savepath=None,debug=0): """ Plot chains to control performance of walkers. """ ndim = len(errors) nsteps = len(samples[:,0,0]) fig, axes = plt.subplots(ndim, figsize=(10, 7), sharex=True) for i in range(ndim): if ndim == 1: ax = axes else: ax = axes[i] ax.fill_between([0.0,nsteps],[-errors[i],-errors[i]],[errors[i],errors[i]],color="g",alpha=0.2) ax.plot(errors[i]*samples[:, :, i], "k", alpha=0.3) ax.axhline(0,c="g",ls=":") ax.set_xlim(0,nsteps) ax.set_ylabel(labels[i]) ax.yaxis.set_label_coords(-0.1, 0.5) if ndim == 1: ax.set_xlabel("step number"); else: axes[-1].set_xlabel("step number"); if savepath is not None: plt.savefig(savepath,bbox_inches='tight') if debug > 1: plt.show() else: plt.show() plt.close() ################################################################################ # FUNCTIONS: Reading files and documentation ################################################################################ def read_input_ft1_file(infile,FT2,weightfield,wmin,planets): "Read FT1 file." try: SatelliteObs(name="Fermi",ft2name=FT2) except ValueError: pass tl = load_Fermi_TOAs(infile,weightcolumn=weightfield,minweight=wmin) toas = toa.get_TOAs_list(tl,planets=planets) weights = np.array([float(w["weight"]) for w in toas.table["flags"]]) energies = np.array([float(e["energy"]) for e in toas.table["flags"]]) return toas, weights, energies def read_pulse_profile_template(prof_file): "Read in LAT pulse profile template." primitives,norms = prim_io(prof_file) template = LCTemplate(primitives,norms) return template def templ_gen(pk_pars,npeaks): "Read template parameters and generate PINT-readable template." tmp_pk_pars = np.copy(pk_pars) pk_phas = tmp_pk_pars[:npeaks] pk_fwhm = tmp_pk_pars[npeaks:-npeaks] pk_ampl = tmp_pk_pars[-npeaks:] if npeaks > 1: pk_phas[1:] = pk_phas[1:] + pk_phas[0] pk_phas = np.mod(pk_phas,1.0) pk_sigm = pk_fwhm / two_sqrt_two_log_two prims = [(LCGaussian(p = [sigma,mu])) for sigma,mu in zip(pk_sigm,pk_phas)] template = LCTemplate(prims,pk_ampl) return template ################################################################################ # FUNCTIONS: Analysis ################################################################################ def check_tau(tau,nburn,nsteps,loud=None): """ [DEPRECATED] Check if some tau values are NaN, and return max(tau). """ if np.isnan(tau).any(): if loud: print("\n# Warning: nan in tau!") print(tau) tmp = np.nanmax(tau) if np.isnan(tmp): if loud: nthin = 1 print("\n# Warning: Actually all values in tau are nan!") else: nthin = np.inf else: nthin = int(np.ceil(tmp)) else: nthin = int(np.ceil(np.max(tau))) if loud: print("\n# Samples thinned by factor {0}".format(nthin)) if float(nburn)/float(nthin) < 10: print("\n# Warning: #burnin < 10 tau.") print(tau) print("# Use more than {0} burnin steps.\n".format(nthin*10)) else: print("# Burnin is {0:.1f} tau".format(float(nburn)/float(nthin))) if float(nsteps)/float(nthin) < 50: print("\n# Warning: #steps < 50 tau.") print(tau) print("# Use more than {0} MCMC steps.\n".format(nthin*60)) else: print("# Sampling is {0:.1f} tau".format(float(nsteps)/float(nthin))) if loud: print("Using tau = {0:d}".format(nthin)) return nthin def mcmc_run(sampler,startpos,nsteps,usetau,doburnin,savepath=None,debug=0): """ Run MCMC analysis calling emcee. """ start = time.time() if doburnin: newpos, newprob, newstate = sampler.run_mcmc(startpos, nsteps, progress=True, store=True) else: if usetau: count = 0 max_n = 100000 goal = 120.0 # (goal x tau) steps tmppos = np.copy(startpos) autocorr,autocorr2 = [],[] alltau,alltau2 = [],[] old_tau = np.inf plt.ion() fig,axes = plt.subplots(nrows=2,ncols=1,sharex=True) for sample in sampler.sample(tmppos, iterations=max_n, progress=True): if sampler.iteration % nsteps: continue tau = np.array(sampler.get_autocorr_time(tol=0)) tau[~np.isfinite(tau)] = 400.0 autocorr.append(np.mean(tau)) autocorr2.append(np.mean(np.ceil(tau))) alltau.append(tau) alltau2.append(np.ceil(tau)) count += 1 # ------ start: show current status if len(autocorr)>1: n2 = nsteps * np.arange(2,count+1) axes[0].clear() for i in range(len(alltau[0])): y2_alltau = np.abs(np.array(alltau)[1:,i]-np.array(alltau)[:-1,i])/np.array(alltau)[:-1,i] axes[0].plot(n2, y2_alltau,c="lightgreen") y2_mean = np.abs(np.array(autocorr)[1:]-np.array(autocorr)[:-1])/np.array(autocorr)[:-1] axes[0].plot(n2,y2_mean,c="g") axes[0].axhline(0.01,ls="--",c="k") axes[0].grid(True) axes[0].set_yscale('log') axes[0].set_ylabel(r"difference $(\hat{\tau}_{n} - \hat{\tau}_{n-1}) / \hat{\tau}_{n}$") n = nsteps * np.arange(1,count+1) axes[1].clear() if len(autocorr)>1: for i in range(len(alltau[0])): axes[1].plot(n, np.array(alltau2)[:,i],c="lightgreen") y_mean = np.array(autocorr2) axes[1].plot(n, y_mean,c="g") axes[1].set_xlim(0, n.max()) if len(autocorr)>10: axes[1].set_ylim(0, 1.05*(np.array(alltau2)[-5:]).max()) else: axes[1].set_ylim(0, 1.05*(np.array(alltau2)).max()) axes[1].plot(n, n / goal, "--k") axes[1].grid(True) axes[1].set_xlabel("number of steps") axes[1].set_ylabel(r"mean $\hat{\tau}$") axes[1].set_title(r"goal: 120$\tau$, current #steps: {0:d}".format(int(nsteps*count))) plt.pause(0.01) plt.draw() plt.pause(0.01) # ------ end: show current status converged = np.all(np.ceil(tau) * goal < sampler.iteration) converged &= np.all(np.abs(old_tau - tau) / tau < 0.01) if converged: break old_tau = tau newpos, newprob, newstate = sampler.get_last_sample() nthin = int(np.ceil(np.max(tau))) axes[0].set_title("converged") plt.ioff() if savepath is not None: plt.savefig(savepath,bbox_inches='tight') if debug > 1: plt.show() else: plt.show() plt.close() else: newpos, newprob, newstate = sampler.run_mcmc(startpos, nsteps, progress=True, store=True) tau = np.array(sampler.get_autocorr_time(tol=0)) tau[~np.isfinite(tau)] = 400.0 nthin = int(np.ceil(np.max(tau))) if nthin == 400: nthin = 1 print("Warning: Not converged!") print("list tau: ", tau) print("discard: ", 20*nthin) print("thin: ", nthin) end = time.time() print("Run took {0:.1f} seconds".format(end-start)) return sampler,newpos,nthin end = time.time() print("Run took {0:.1f} seconds".format(end-start)) return sampler,newpos ################################################################################ # FUNCTIONS: Test statistics ################################################################################ def likelihood(phases,weights,template): """ Calculate likelihood at all phase shifts between 0 and 1, summed over all photons in segment. """ like = np.prod((1.0-weights)+weights*template(phases)) return like def loglike(phases,weights,template,dphi=0.0): """ Calculate log-likelihood at all phase shifts between 0 and 1, summed over all photons in segment. """ llike = np.sum(np.log((1.0-weights)+weights*template(np.mod(phases-dphi,1.0)))) return llike def htest(phases,weights): """ H statistic for given phases and weights. """ cosSum = np.cumsum(weights * np.cos(2.0 * np.pi * np.arange(1,21)[np.newaxis].T * phases),axis=1) sinSum = np.cumsum(weights * np.sin(2.0 * np.pi * np.arange(1,21)[np.newaxis].T * phases),axis=1) norm = 2.0 / np.cumsum(weights ** 2.0) PP = ((cosSum ** 2.0) + (sinSum ** 2.0)) * norm Htest = np.max(np.cumsum(PP,axis=0) - 4.0 * np.arange(1,21)[np.newaxis].T + 4.0, axis=0) return Htest[-1] def chi_square(res,sig): """ Chi-squared statistic. """ mean = np.average(res,weights=1/(sig*sig)) return np.sum(((res-mean)/sig)**2) def red_chi_square(res,sig,n=0): """ Reduced chi-squared statistic. """ nres = len(res) return chi_square(res,sig)/(nres-n) ################################################################################ # FUNCTIONS: Joint radio and gamma-ray test statistics ################################################################################ def log_prior(values,theta,labels): """ Priors set to -inf outside of given parameter range. Other priors possible. """ # Sanity checks for parameters if len(theta)>=1 and theta[0]=="OFFSET": if -0.5 > theta[0] >= 0.5: return -np.inf if "F1" in labels: id = np.where(labels=="F1")[0][0] if values[id]+theta[id]>0.0: return -np.inf if "A1" in labels: id = np.where(labels=="A1")[0][0] if values[id]+theta[id]<0.0: return -np.inf return 0.0 def log_prob_precise(theta,fitpars,labels,values,ranges,npeaks,gt,gw,rt,rre,tm_base): """ Compute joint log-likelihood using PINT. Input: -- Theta, gamma-toas, gamma-weights, radio-toas, radio-residual-errors, timing model -- """ # multiply errors given in parfile with theta paroffsets = theta * ranges if npeaks is not None: # split input theta in pulsar pars and template profile pars labs = labels[:-npeaks*3] pars = values[:-npeaks*3] + paroffsets[:-npeaks*3] temp = values[-npeaks*3:] + paroffsets[-npeaks*3:] if len(np.where((temp[npeaks:-npeaks] < lob_fwhm) | (temp[npeaks:-npeaks] > hib_fwhm))[0]) > 0: return -np.inf if len(np.where(temp[-npeaks:] <= 0.0)[0]) > 0: return -np.inf if np.sum(temp[-npeaks:]) > 1: return -np.inf tmp_template = templ_gen(temp,npeaks) else: labs = labels pars = values + paroffsets # read in pulsar pars tm = tm_base if labels[0] != "OFFSET": for i,p in enumerate(fitpars): getattr(tm,p).value = pars[i] offset = 0.0 else: for i,p in enumerate(fitpars): getattr(tm,p).value = pars[i+1] offset = pars[0] # set prior for loglike lp = log_prior(values,pars,labs) if not np.isfinite(lp): # go home if prior doesnt like theta return -np.inf # compute llike for gamma data if gt is None: # no gamma data -> L = 0 gllike = 0.0 else: if npeaks is None: # w/o template no real llike gp = np.mod(tm.phase(gt).frac,1.0).value gllike = 0.4*htest(gp,gw) else: gp = np.mod(tm.phase(gt).frac,1.0) gllike = loglike(gp,gw,tmp_template)#,dphi=offset) # compute llike for radio data if rt is None: # no radio data -> L = 0 rllike = 0.0 else: rr = Residuals(rt,tm).phase_resids rllike = -0.5*chi_square(rr,rre) return lp+gllike+rllike def log_prob_quick(theta,labels,values,ranges,npeaks,gM,gp,gw,rM,rr,rre): """ Compute joint log-likelihood using PINT. Input: -- Theta, gamma-phases, gamma-weights, radio-residuals, radio-residual-error, model-derivatives -- """ # multiply errors given in parfile with theta paroffsets = theta * ranges if npeaks is not None: # split input theta in pulsar pars and template profile pars labs = labels[:-npeaks*3] pars = paroffsets[:-npeaks*3] temp = values[-npeaks*3:] + paroffsets[-npeaks*3:] if len(np.where((temp[npeaks:-npeaks] < lob_fwhm) | (temp[npeaks:-npeaks] > hib_fwhm))[0]) > 0: #print(temp[npeaks:-npeaks]) #print("FWHM") return -np.inf if len(np.where(temp[-npeaks:] <= 0.0)[0]) > 0: #print("AMP 0") return -np.inf if np.sum(temp[-npeaks:]) > 1: #print("AMP 1") return -np.inf tmp_template = templ_gen(temp,npeaks) else: if labels[0] == "OFFSET": labs = labels[1:] pars = paroffsets[1:] else: labs = labels pars = paroffsets # set prior for loglike lp = log_prior(values,pars,labs) if not np.isfinite(lp): # go home if prior doesnt like theta #print("PRIOR") return -np.inf # compute llike for gamma data if gp is None: # no gamma data -> L = 0 gllike = 0.0 else: gdphi = np.matmul(gM,pars) if npeaks is None: # w/o template no real llike gllike = 0.4*htest(gp-gdphi,gw) else: gllike = loglike(gp-gdphi,gw,tmp_template) # compute llike for radio data if rr is None: # no radio data -> L = 0 rllike = 0.0 else: rdphi = np.matmul(rM,pars) rllike = -0.5*chi_square(rr-rdphi,rre) return lp+gllike+rllike ################################################################################ # MAIN: ################################################################################ if __name__ == '__main__': ############################################################################ # INPUT: Command line options ############################################################################ description = '##########################################################\ # Joint radio and gamma-ray pulsar timing (jrag timing). #\ ##########################################################' parser = OptionParser(usage = ' %prog [options]', description = description) baseparser = OptionGroup(parser, "Basic input parameters.", "Caution: Some of these are required.") baseparser.add_option('-p','--parfile',type='str',default=None,help='Path to parameter file [REQUIRED].') baseparser.add_option('--newparfile',type='str',default=None,help='Output parameter file.') baseparser.add_option('--resultsfile',type='str',default=None,help='Output results file.') baseparser.add_option('-t','--tempfile',type='str',default=None,help='Path to template file.') baseparser.add_option('--newtempfile',type='str',default=None,help='Output template file.') baseparser.add_option('-l','--ft1file',type='str',default=None,help='Path to FT1 file (LAT data).') baseparser.add_option('-f','--ft2file',type='str',default=None,help='Path to FT2 file (satellite).') baseparser.add_option('-r','--timfile',type='str',default=None,help='Path to tim file (radio).') baseparser.add_option('--weightfield',type='str',default='Weights',help='Weightfield name in FT1 file.') baseparser.add_option('-d','--debug',type='int',default=1,help='Level of debug messages.') baseparser.add_option('--plotonly',type='int',default=0,help='No MCMC run, only plotting from previous runs.') parser.add_option_group(baseparser) timparser = OptionGroup(parser, "Timing options", "Caution: Choose wisely!") # total number of walkers = (multiple of 2x number of CPU threads) timparser.add_option('-n','--npoints',type='int',default=128,help='#walkers. [default: %default]') # maximum number of iterations - goal: 100x autocorr time (!! overwritten by --usetau=1) timparser.add_option('-x','--maxiter',type='int',default=100,help='#iterations. [default: %default]') # number of burnin iterations - goal: 20x autocorr time (!! gets redone if --usetau=1) timparser.add_option('-b','--burnin',type='int',default=100,help='#burnins. [default: %default]') # use autocorr time (tau) for #burnin and #nsteps timparser.add_option('--usetau',type='int',default=0,help='#tau in timing. [default: %default]') # initial offset of parameters from previous best solution timparser.add_option('-s','--offset',type='float',default=0.01,help='Fractional offset. [default: %default]') # photon weight cutoff timparser.add_option('-c','--pcut',type='float',default=0.0,help='Minimum weight to be used. [default: %default]') # photon weight cutoff timparser.add_option('--eclipse',type='float',default=0.0,help='Orbital phase radius in which radio TOAs are removed. [default: %default]') # photon weight cutoff timparser.add_option('--toacut',type='float',default=1000.0,help='Maximum error in radio TOAs (microsecs). [default: %default]') # compare radio and gamma-ray posteriors timparser.add_option('--compare',type='int',default=0,help='Compare radio and gamma-ray posteriors. [default: %default]') parser.add_option_group(timparser) (options,args) = parser.parse_args() if options.debug > 2: print("Options:") print(options) print("Args:") print(args) ############################################################################ # INPUT: warning and error messages for input ############################################################################ if options.parfile == None: print("### Error: No parfile specified!", file=sys.stderr) print("Use -p Jxxxx-xxxx_n1_v1.par to specify the parfile.", file=sys.stderr) sys.exit(1) else: if not os.path.isfile(options.parfile): print("### Error: Specified parfile does not exist!", file=sys.stderr) sys.exit(1) if ".pa" in options.parfile and "_v" in options.parfile and "_n" in options.parfile: base, number = options.parfile.rsplit("_v", 1)[0].rsplit("_n", 1) tmp, end = options.parfile.split("_n{0:s}".format(number)) source = base.rsplit("/", 1)[1].split("_", 1)[0] if options.newparfile == None: if ".pa" in options.parfile and "_v" in options.parfile: base, number = options.parfile.rsplit(".pa", 1)[0].rsplit("_v", 1) options.newparfile = "{0:s}_v{1:d}.par".format(base, int(number) + 1) else: print("No name for new parfile specified!", file=sys.stderr) print("Use --newparfile Jxxxx-xxxx_n1_v2.par to specify the new parfile.", file=sys.stderr) sys.exit(1) if options.resultsfile == None: if ".pa" in options.parfile and "_v" in options.parfile: base, number = options.parfile.rsplit(".pa", 1)[0].rsplit("_v", 1) options.resultsfile = "{0:s}_v{1:d}.dat".format(base, int(number)) else: print("No name for resultsfile specified!", file=sys.stderr) print("Use --resultsfile Jxxxx-xxxx_n1_v1.dat to specify the resultsfile.", file=sys.stderr) sys.exit(1) if options.ft1file == None or options.ft2file == None: print("## Warning: No FT1 and/or FT2 file given. Will do no gamma-ray timing!") dogamma = 0 else: if not os.path.isfile(options.ft1file): print("### Error: Specified FT1 file does not exist!", file=sys.stderr) sys.exit(1) if not os.path.isfile(options.ft2file): print("### Error: Specified FT2 file does not exist!", file=sys.stderr) sys.exit(1) dogamma = 1 # only for gamma timing a template file is needed if options.tempfile == None: if ".pa" in options.parfile and "_v" in options.parfile: base, number = options.parfile.rsplit(".pa", 1)[0].rsplit("_v", 1) options.tempfile = "{0:s}_v{1:d}.itemp".format(base, int(number)) if not os.path.isfile(options.tempfile): print("### Error: Expected template file does not exist!", file=sys.stderr) sys.exit(1) else: print("No name for tempfile specified!", file=sys.stderr) print("Use --tempfile Jxxxx-xxxx_n1_v1.itemp to specify the resultsfile.", file=sys.stderr) sys.exit(1) else: if not os.path.isfile(options.tempfile): print("### Error: Specified template file does not exist!", file=sys.stderr) sys.exit(1) # see above if options.newtempfile == None: if ".pa" in options.parfile and "_v" in options.parfile: base, number = options.parfile.rsplit(".pa", 1)[0].rsplit("_v", 1) options.newtempfile = "{0:s}_v{1:d}.itemp".format(base, int(number) + 1) else: print("No name for new parfile specified!", file=sys.stderr) print("Use --newparfile Jxxxx-xxxx_n1_v2.itemp to specify the new parfile.", file=sys.stderr) sys.exit(1) if options.timfile == None: print("## Warning: No timfile given. Will do no radio timing!") doradio = 0 else: if not os.path.isfile(options.timfile): print("### Error: Specified timfile does not exist!", file=sys.stderr) sys.exit(1) doradio = 1 if dogamma+doradio == 2: doboth = 1 else: doboth = 0 if options.compare == 1: print("### Error: Compare = 1, but not both data given.", file=sys.stderr) sys.exit(1) if options.plotonly == 0: # Start writing in log file logfile = open("{0:s}.log".format(options.resultsfile),'w') logfile.write("#-------------------------------------------------------------------------------------------------\n") logfile.write("# Logfile\n") logfile.write("#-------------------------------------------------------------------------------------------------\n\n") line = "python" for input_str in sys.argv: line += " {0:s}".format(input_str) logfile.write(line+"\n\n") ############################################################################ # PREPARATION: Read in timing model and template file ############################################################################ timing_model = get_model(options.parfile) if options.debug > 0: print("\n# Timing model:") print(timing_model) if getattr(timing_model,"PLANET_SHAPIRO").quantity == True: planets = True else: planets = False if getattr(timing_model,"BINARY").quantity == None: isolated = True else: isolated = False fitpars = np.array([p for p in timing_model.params if not getattr(timing_model,p).frozen]) # sort out all jumps jumppars,jumpvals,jumperrs = [],[],[] for par_name in fitpars: if "JUMP" in par_name: jumppars.append(par_name) if options.debug > 3: print(getattr(timing_model,par_name)) #print(getattr(timing_model,par_name).key_value) jumpunits = ["s"]*len(jumppars) if len(jumppars) > 0: jumpvals = [getattr(timing_model,p).value for p in jumppars] jumperrs = [getattr(timing_model,p).uncertainty_value for p in jumppars] jumptels = [getattr(timing_model,p).key_value[0] for p in jumppars] if options.debug > 1: print("JUMP pars:") print(jumppars) print(jumpvals) print(jumperrs) print(jumptels) print(jumpunits) # put parameters in known order possible_ordered_pars = np.array(["RAJ","DECJ","PMRA","PMDEC","PX","F0","F1","F2","F3","F4","F5","DM","PB","A1","TASC","EPS1","EPS2","ECC","OM","PBDOT","FB1","FB2","FB3","FB4","FB5"]+jumppars) possible_ordered_units = np.array(["rad","rad","mas/yr","mas/yr","mas","Hz","Hz/s","Hz/s^2","Hz/s^3","Hz/s^4","Hz/s^5","1","d","lts","d","1","1","1","deg","d/d","Hz/s","Hz/s^2","Hz/s^3","Hz/s^4","Hz/s^5"]+jumpunits) if options.debug > 0: print("\n# Pulsar parameters to fit:") print(fitpars) order_mask = [] fitunits = [] for par_name in possible_ordered_pars: if par_name in fitpars: order_mask.append(np.where(fitpars==par_name)[0][0]) fitunits.append(possible_ordered_units[np.where(possible_ordered_pars==par_name)[0][0]]) fitpars = (fitpars[order_mask]).tolist() if options.debug > 0: print("# Sorted (w/o jumps).") print(fitpars) fitvals = [getattr(timing_model,p).value for p in fitpars] fiterrs = [getattr(timing_model,p).uncertainty_value for p in fitpars] if options.debug > 1: print(fitvals) print(fiterrs) npeaks = None if options.tempfile != None: npeaks = 0 # latprofile = read_pulse_profile_template(options.tempfile) with open(options.tempfile,"r") as tfile: lines = tfile.readlines() #print(lines) npeaks_check = 0 for line in lines: if "phas" in line: npeaks += 1 if int(line.split()[0][-1]) > npeaks_check: npeaks_check = int(line.split()[0][-1]) if npeaks != npeaks_check: print("### Error: Number of peaks and parameter names are inconsistent.", file=sys.stderr) print("Number of peaks: {0}".format(npeaks)) print("Highest number in peak name: {0}".format(npeaks_check)) sys.exit(1) pk_phas_val = np.zeros(npeaks) pk_phas_err = np.zeros(npeaks) pk_fwhm_val = np.zeros(npeaks) pk_fwhm_err = np.zeros(npeaks) pk_ampl_val = np.zeros(npeaks) pk_ampl_err = np.zeros(npeaks) for line in lines: splittedline = line.strip().split() if "phas" in line: idx = int(splittedline[0][-1]) pk_phas_val[idx-1] = float(splittedline[2]) pk_phas_err[idx-1] = float(splittedline[4]) elif "fwhm" in line: idx = int(splittedline[0][-1]) pk_fwhm_val[idx-1] = float(splittedline[2]) pk_fwhm_err[idx-1] = float(splittedline[4]) elif "ampl" in line: idx = int(splittedline[0][-1]) pk_ampl_val[idx-1] = float(splittedline[2]) pk_ampl_err[idx-1] = float(splittedline[4]) if npeaks > 1: if options.debug > 2: print("Profile peak phases:") print(pk_phas_val) print(pk_phas_val[1:]) print(pk_phas_val[1:] - pk_phas_val[0]) print("Profile peak amplitudes:") print(pk_ampl_val) pk_phas_val[1:] = pk_phas_val[1:] - pk_phas_val[0] pk_pars = np.concatenate((pk_phas_val,pk_fwhm_val,pk_ampl_val)) latprofile_orig = templ_gen(pk_pars,npeaks) if npeaks == 1: pk_labels = ["PHI0","FWHM0","AMPL0"] elif npeaks > 1: pk_labels = ["PHI0"] for i in range(1,npeaks): pk_labels.append("DPHI{0}".format(i)) for i in range(npeaks): pk_labels.append("FWHM{0}".format(i)) for i in range(npeaks): pk_labels.append("AMPL{0}".format(i)) print(pk_labels) pk_units = ["1"]*3*npeaks if options.debug > 0 or not latprofile_orig.norm_ok(): print("\n# Input LAT pulse profile template:") print(latprofile_orig) print("Norm: {0}".format(latprofile_orig.norm())) # arrays required by logL-functions if doboth: #labels = np.concatenate((["OFFSET"],fitpars,pk_labels)) #offsetgate #values = np.concatenate(([0.0],fitvals,pk_pars)) #errors = np.concatenate(([0.01],fiterrs,pk_phas_err,pk_fwhm_err,pk_ampl_err)) #units = np.concatenate((["1"],fitunits,pk_units)) labels = np.concatenate((fitpars,pk_labels)) values = np.concatenate((fitvals,pk_pars)) errors = np.concatenate((fiterrs,pk_phas_err,pk_fwhm_err,pk_ampl_err)) junits = np.concatenate((fitunits,pk_units)) elif doradio: labels = np.copy(fitpars) values = np.copy(fitvals) errors = np.copy(fiterrs) junits = np.copy(fitunits) elif dogamma: if len(jumppars) > 0: labels = np.concatenate((fitpars[:-len(jumppars)],pk_labels)) values = np.concatenate((fitvals[:-len(jumppars)],pk_pars)) errors = np.concatenate((fiterrs[:-len(jumppars)],pk_phas_err,pk_fwhm_err,pk_ampl_err)) junits = np.concatenate((fitunits[:-len(jumppars)],pk_units)) else: labels = np.concatenate((fitpars,pk_labels)) values = np.concatenate((fitvals,pk_pars)) errors = np.concatenate((fiterrs,pk_phas_err,pk_fwhm_err,pk_ampl_err)) junits = np.concatenate((fitunits,pk_units)) nfitpars = len(labels) if options.debug > 1: print("Number of fit pars: {0:d}".format(nfitpars)) for i in range(nfitpars): print("{0:<7} {1:<8.3g} {2:.3g} {3:s}".format(labels[i],values[i],errors[i],junits[i])) if doradio: possible_radio_pars = np.copy(possible_ordered_pars) mask_radio_pars = [] for i in range(len(labels)): if labels[i] in possible_radio_pars: mask_radio_pars.append(i) rlabels = labels[mask_radio_pars] if dogamma: mask = [] radio_only_pars = np.array(["DM"]+jumppars) for i in range(len(possible_ordered_pars)): if possible_ordered_pars[i] not in radio_only_pars: mask.append(i) possible_gamma_pars = np.concatenate((possible_ordered_pars[mask],pk_labels)) mask_gamma_pars = [] for i in range(len(labels)): if labels[i] in possible_gamma_pars: mask_gamma_pars.append(i) glabels = labels[mask_gamma_pars] if options.compare: compareskylabels = np.array(["RAJ","DECJ","PMRA","PMDEC","F0","F1","F2","PB","A1","TASC","EPS1","EPS2"]) mask_compare_pars = [] for i in range(len(compareskylabels)): if compareskylabels[i] in labels: mask_compare_pars.append(i) compareskylabels = np.copy(compareskylabels[mask_compare_pars]) if len(compareskylabels) == 0: print("Attention: No compare pars.") mask_compare_both = [] for i in range(len(labels)): if labels[i] in compareskylabels: mask_compare_both.append(i) mask_compare_radio = [] for i in range(len(rlabels)): if rlabels[i] in compareskylabels: mask_compare_radio.append(i) mask_compare_gamma = [] for i in range(len(glabels)): if glabels[i] in compareskylabels: mask_compare_gamma.append(i) if options.debug > 1: if doradio: print("Number of radio fit pars: {0:d}".format(len(mask_radio_pars))) print(rlabels) if dogamma: print("Number of gamma fit pars: {0:d}".format(len(mask_gamma_pars))) print(glabels) ############################################################################ # PREPARATION: Read in data and calculate test statistics ############################################################################ # --------------------------------------------------- # Radio TOAs if doradio: radio_toas_orig = toa.get_TOAs(options.timfile,planets=planets) radio_residuals = Residuals(radio_toas_orig,timing_model).phase_resids if options.debug > 1: print("Original radio TOAs:\n") radio_toas_orig.print_summary() plot_radio_residuals(radio_toas_orig,timing_model,isolated=isolated,savepath="{0:s}_cleanRes1.png".format(options.resultsfile),debug=options.debug) if not options.toacut and not options.eclipse: if options.debug > 1: print("No TOA cleaning done. Default values:") print("[--toacut 1000.0 (error in microseconds)]") print("[--eclipse 0.0 (eclipsed orb. phase)]") radio_toas = radio_toas_orig else: # clean TOAs from large uncertainties error_range = options.toacut * u.us error_notok = radio_toas_orig.table["error"] > error_range orbphi_range = options.eclipse if isolated == False: # clean TOAs from orbital phases ... # potentially affected by intrabinary material radio_orbphi = timing_model.orbital_phase(radio_toas_orig,anom="mean",radians=False) orbphi_notok = (radio_orbphi > 0.25-0.5*orbphi_range) & (radio_orbphi < 0.25+0.5*orbphi_range) else: orbphi_notok = np.array([False]*len(error_notok)) # final cleaned TOAs radio_toas = radio_toas_orig[~(error_notok+orbphi_notok)] if options.debug > 0: print("Cleaned TOAs:\n") radio_toas.print_summary() if options.debug > 1: if (error_notok+orbphi_notok).sum() == 0: print("No TOAs removed.\n") else: print("Removed TOAs:\n") (radio_toas_orig[error_notok+orbphi_notok]).print_summary() plot_radio_residuals(radio_toas_orig,timing_model,orbphi_range=orbphi_range,error_range=error_range,isolated=isolated,savepath="{0:s}_cleanRes2.png".format(options.resultsfile)) plot_radio_residuals(radio_toas,timing_model,orbphi_range=orbphi_range,isolated=isolated,savepath="{0:s}_cleanRes.png".format(options.resultsfile),debug=options.debug) radio_toa_uncerts = (radio_toas.get_errors()*timing_model.F0.quantity).to_value(u.dimensionless_unscaled) radio_residuals = Residuals(radio_toas,timing_model).phase_resids M_radio,names,units = timing_model.designmatrix(radio_toas) if doboth: M_radio = M_radio[:,(np.array(order_mask)+1).tolist()] else: M_radio = M_radio[:,(np.array(order_mask)+1).tolist()] M_radio *= timing_model.F0.value else: radio_toas = None radio_toa_uncerts = None radio_residuals = None M_radio = None # --------------------------------------------------- # Gamma rays if dogamma: fermi_toas, weights, energies = read_input_ft1_file(options.ft1file,options.ft2file,options.weightfield,options.pcut,planets) phases = np.mod(timing_model.phase(fermi_toas).frac,1.0) M_gamma,names,units = timing_model.designmatrix(fermi_toas) if doboth: M_gamma = M_gamma[:,(np.array(order_mask)+1).tolist()] else: if len(jumppars) > 0: M_gamma = M_gamma[:,(np.array(order_mask)+1).tolist()][:,:-len(jumppars)] else: M_gamma = M_gamma[:,(np.array(order_mask)+1).tolist()] M_gamma *= timing_model.F0.value if options.debug > 0: print("Scanning phase for PHI0...") if errors[-npeaks*3] == 0.0: print("WARNING: No uncertainty given on template gamma profile.") reldphi_list = np.linspace(-0.5/errors[-npeaks*3],0.5/errors[-npeaks*3],200) l2dphi_list,l2gdphi_list,l2rdphi_list = [],[],[] for i in range(len(reldphi_list)): theta = np.zeros(nfitpars) theta[-npeaks*3] = reldphi_list[i] l2dphi_list.append(log_prob_quick(theta,labels,values,errors,npeaks,M_gamma,phases,weights,M_radio,radio_residuals,radio_toa_uncerts)) l2gdphi_list.append(log_prob_quick(theta,labels,values,errors,npeaks,M_gamma,phases,weights,None,None,None)) l2rdphi_list.append(log_prob_quick(theta,labels,values,errors,npeaks,None,None,None,M_radio,radio_residuals,radio_toa_uncerts)) if options.debug > 2: fig, ax = plt.subplots(3, figsize=(10, 7), sharex=True) ax[0].plot(reldphi_list*errors[-npeaks*3],l2dphi_list,c="k",ls="-") ax[1].plot(reldphi_list*errors[-npeaks*3],l2gdphi_list,c="b",ls="-.") ax[2].plot(reldphi_list*errors[-npeaks*3],l2rdphi_list,c="r",ls="--") plt.savefig("{0:s}_PHI0scan.png".format(options.resultsfile),bbox_inches='tight') plt.show() plt.close() plot_phases(fermi_toas,weights,radio_toas,timing_model,isolated=isolated,template=latprofile_orig) pk_pars[0] += errors[-npeaks*3]*reldphi_list[np.where(l2gdphi_list==np.max(l2gdphi_list))][0] values[-npeaks*3] += errors[-npeaks*3]*reldphi_list[np.where(l2gdphi_list==np.max(l2gdphi_list))][0] latprofile = templ_gen(pk_pars,npeaks) if options.debug > 2: plot_phases(fermi_toas,weights,radio_toas,timing_model,isolated=isolated,template=latprofile) if options.debug > 0: print("Rotated template: {0:.4f}".format(errors[-npeaks*3]*reldphi_list[np.where(l2gdphi_list==np.max(l2gdphi_list))][0])) print("Scanning phase for PHI0... Done.") else: fermi_toas = None phases = None weights = None M_gamma = None if options.debug > 0: theta = np.zeros(nfitpars) print("\n---------------------------------------------------------------\n") print("Test statistics for testing:") if doradio: print("\nChiSquared: {0: 14.2f}".format(chi_square(radio_residuals,radio_toa_uncerts))) print("LLike from ChiSquared: {0: 14.2f}".format(-0.5*chi_square(radio_residuals,radio_toa_uncerts))) if dogamma: print("\nLogLikelihood: {0: 14.2f}".format(loglike(phases,weights,latprofile))) print("H-test: {0: 14.2f}".format(htest(phases.value,weights))) print("LLike from H-test: {0: 14.2f}".format(0.4*htest(phases.value,weights))) if doboth: print("\nJoint Radio and Gamma LogL using H-test or LogL (for gamma) and Chi2 (for radio):") print("Joint H test: {0: 14.2f}".format(0.4*htest(phases.value,weights)\ -0.5*chi_square(radio_residuals,radio_toa_uncerts))) #print("Using function: {0: 14.2f}".format(log_prob_precise(theta,fitpars,values,errors,npeaks,fermi_toas,weights,radio_toas,radio_toa_uncerts,timing_model))) print("Joint LogLike: {0: 14.2f}".format(loglike(phases,weights,latprofile)\ -0.5*chi_square(radio_residuals,radio_toa_uncerts))) print("Using function: {0: 14.2f}".format(log_prob_precise(theta,fitpars,labels,values,errors,npeaks,fermi_toas,weights,radio_toas,radio_toa_uncerts,timing_model))) print("Using fast function: {0: 14.2f}".format(log_prob_quick(theta,labels,values,errors,npeaks,M_gamma,phases,weights,M_radio,radio_residuals,radio_toa_uncerts))) print("\n---------------------------------------------------------------\n") if options.debug > 3 and "F1" in labels and "DECJ" in labels and options.plotonly==0: # prepare search over F1 nfdot = 100 id_fdot = np.where(labels == "F1")[0][0] dfdot_list = np.linspace(-1.0,1.0,num=nfdot)*3 fdot_step = (dfdot_list[1]-dfdot_list[0])*errors[id_fdot] # prepare search over DECJ ndecj = 100 id_decj = np.where(labels == "DECJ")[0][0] ddecj_list = np.linspace(-1.0,1.0,num=ndecj)*3 decj_step = (ddecj_list[1]-ddecj_list[0])*errors[id_decj] ljgrid_list = np.zeros((len(dfdot_list),len(ddecj_list))) lggrid_list = np.zeros((len(dfdot_list),len(ddecj_list))) lrgrid_list = np.zeros((len(dfdot_list),len(ddecj_list))) for i in range(len(ddecj_list)): theta = np.zeros(nfitpars) theta[id_decj] = ddecj_list[i] for ii in range(len(dfdot_list)): theta[id_fdot] = dfdot_list[ii] ljgrid_list[ii,i] = log_prob_quick(theta,labels,values,errors,npeaks,M_gamma,phases,weights,M_radio,radio_residuals,radio_toa_uncerts) lggrid_list[ii,i] = log_prob_quick(theta,labels,values,errors,npeaks,M_gamma,phases,weights,None,None,None) lrgrid_list[ii,i] = log_prob_quick(theta,labels,values,errors,npeaks,None,None,None,M_radio,radio_residuals,radio_toa_uncerts) scan_Lmax = np.max(ljgrid_list) scan_max = np.where(ljgrid_list == scan_Lmax) extent = [np.min(ddecj_list)*errors[id_decj]-0.5*decj_step,np.max(ddecj_list)*errors[id_decj]+0.5*decj_step,\ np.min(dfdot_list)*errors[id_fdot]-0.5*fdot_step,np.max(dfdot_list)*errors[id_fdot]+0.5*fdot_step] ddecj_max = ddecj_list[scan_max[1][0]]*errors[id_decj] dfdot_max = dfdot_list[nfdot-1-scan_max[0][0]]*errors[id_fdot] decj_max = values[id_decj] + ddecj_max fdot_max = values[id_fdot] + dfdot_max fig, ax = plt.subplots(3, figsize=(10, 7), sharex=True, sharey=True) Lrange = 12.5 c0=ax[0].imshow(ljgrid_list,extent=extent,aspect="auto",vmin=scan_Lmax-Lrange) ax[0].scatter(ddecj_max,dfdot_max,marker="+",color="r") ax[0].set_ylabel("Delta F1") cb0 = plt.colorbar(c0) cb0.ax.set_ylabel("L_joint") c1=ax[1].imshow(lggrid_list,extent=extent,aspect="auto",vmin=np.max(lggrid_list)-Lrange) ax[1].scatter(ddecj_max,dfdot_max,marker="+",color="r") ax[1].set_ylabel("Delta F1") cb1 = plt.colorbar(c1) cb1.ax.set_ylabel("L_gamma") c2=ax[2].imshow(lrgrid_list,extent=extent,aspect="auto",vmin=np.max(lrgrid_list)-Lrange) ax[2].scatter(ddecj_max,dfdot_max,marker="+",color="r") ax[2].set_ylabel("Delta F1") ax[2].set_xlabel("Delta DECJ") cb2 = plt.colorbar(c2) cb2.ax.set_ylabel("L_radio") plt.savefig("{0:s}_F1DECJscan.png".format(options.resultsfile),bbox_inches='tight') plt.show() plt.close() # print out results print("\nResults of 2D-scan in F1 and DECJ:") print("Max joint logL: {0:.2f}".format(scan_Lmax)) decj_tmp = Angle(decj_max, u.degree) if decj_tmp.dms.d < 0.0: parvalue = "-" else: parvalue = "+" parvalue += "{0:02.0f}:".format(np.abs(decj_tmp.dms.d)) parvalue += "{0:02.0f}:".format(np.abs(decj_tmp.dms.m)) if np.abs(decj_tmp.dms.s) < 10.0: parvalue += "0{0:.15f}".format(np.abs(decj_tmp.dms.s)) else: parvalue += "{0:.15f}".format(np.abs(decj_tmp.dms.s)) print("Max logL DECJ: {0:s}".format(parvalue)) parvalue = "{0:.16g}".format(fdot_max) print("Max logL F1: {0:s}\n".format(parvalue)) ############################################################################ # MAIN: MCMC run ############################################################################ nwalker = options.npoints nburn = options.burnin if options.maxiter: nsteps = options.maxiter else: nsteps = 100 pos = (options.offset*np.random.randn(nwalker-1,nfitpars)).tolist() pos_null = (np.zeros(nfitpars)).tolist() pos.append(pos_null) # include current best fit pos = np.array(pos) #------------------------------------------------------------ if doboth: if options.plotonly==0: startpos = np.copy(pos) sampler = emcee.EnsembleSampler( nwalker,nfitpars,log_prob_quick, args=(labels,values,errors,npeaks,M_gamma,phases,weights,M_radio,radio_residuals,radio_toa_uncerts), moves=[(emcee.moves.KDEMove(), 1.0),],) if nburn > 0: sampler,burnpos = mcmc_run(sampler,startpos,nburn,usetau=0,doburnin=1) samples = sampler.get_chain() plot_chains(samples,errors,labels, savepath="{0:s}_{1:s}.png".format(options.resultsfile,"burnin_jrag"),debug=options.debug) else: burnpos = np.copy(startpos) sampler.reset() sampler,jragpos,nthin = mcmc_run(sampler,burnpos,nsteps,usetau=options.usetau,doburnin=0, savepath="{0:s}_{1:s}.png".format(options.resultsfile,"converg_jrag"),debug=options.debug) samples = sampler.get_chain() if options.usetau: flat_samples = sampler.get_chain(discard=20*nthin, thin=nthin, flat=True) log_prob_samples = sampler.get_log_prob(discard=20*nthin, thin=nthin, flat=True) else: flat_samples = sampler.get_chain(discard=0, thin=nthin, flat=True) log_prob_samples = sampler.get_log_prob(discard=0, thin=nthin, flat=True) # store samples in compressed file samplesdict = {} samplesdict["samples"] = samples samplesdict["flat_samples"] = flat_samples samplesdict["log_prob_samples"] = log_prob_samples np.savez_compressed("{0:s}_chain".format(options.resultsfile), **samplesdict) loaded = np.load("{0:s}_chain.npz".format(options.resultsfile)) samples = loaded["samples"] flat_samples = loaded["flat_samples"] log_prob_samples = loaded["log_prob_samples"] plot_chains(samples,errors,labels, savepath="{0:s}_{1:s}.png".format(options.resultsfile,"samples_jrag"),debug=options.debug) if options.debug > 3: print("\n # Chain dictionary, loaded keys") print(samplesdict) print(list(loaded.keys())) if options.debug > 2: print("\n# Shape of samples, shape of L samples, maxL, maxL-values") print(flat_samples.shape) print(log_prob_samples.shape) print(np.max(log_prob_samples)) print(flat_samples[np.where(log_prob_samples == np.max(log_prob_samples))][0]) if doradio+doboth == 1 or options.compare == 1: radiolabels = labels[mask_radio_pars] if options.plotonly==0: if doboth == 0: startpos = np.copy(pos) else: startpos = np.copy(jragpos[:,mask_radio_pars]) radiosampler = emcee.EnsembleSampler( nwalker,len(mask_radio_pars),log_prob_quick, #args=(labels, values, errors, npeaks,M_gamma,phases,weights,M_radio, radio_residuals,radio_toa_uncerts), args =(radiolabels,values[mask_radio_pars],errors[mask_radio_pars],None ,None, None, None, M_radio[:,mask_radio_pars],radio_residuals,radio_toa_uncerts), moves=[(emcee.moves.KDEMove(), 1.0),],) if nburn > 0: radiosampler,rburnpos = mcmc_run(radiosampler,startpos,nburn,usetau=0,doburnin=1) radiosamples = radiosampler.get_chain() plot_chains(radiosamples, errors[mask_radio_pars],labels[mask_radio_pars], savepath="{0:s}_{1:s}.png".format(options.resultsfile,"burnin_radio"),debug=options.debug) else: rburnpos = np.copy(startpos) radiosampler.reset() radiosampler,radiopos,nrthin = mcmc_run(radiosampler,rburnpos,nsteps,usetau=options.usetau,doburnin=0, savepath="{0:s}_{1:s}.png".format(options.resultsfile,"converg_radio"),debug=options.debug) radiosamples = radiosampler.get_chain() if options.usetau: flat_radiosamples = radiosampler.get_chain(discard=20*nrthin, thin=nrthin, flat=True) log_prob_radiosamples = radiosampler.get_log_prob(discard=20*nrthin, thin=nrthin, flat=True) else: flat_radiosamples = radiosampler.get_chain(discard=0, thin=nrthin, flat=True) log_prob_radiosamples = radiosampler.get_log_prob(discard=0, thin=nrthin, flat=True) # store samples in compressed file radiosamplesdict = {} radiosamplesdict["radiosamples"] = radiosamples radiosamplesdict["flat_radiosamples"] = flat_radiosamples radiosamplesdict["log_prob_radiosamples"] = log_prob_radiosamples np.savez_compressed("{0:s}_chain_radio".format(options.resultsfile), **radiosamplesdict) loaded = np.load("{0:s}_chain_radio.npz".format(options.resultsfile)) radiosamples = loaded["radiosamples"] flat_radiosamples = loaded["flat_radiosamples"] log_prob_radiosamples = loaded["log_prob_radiosamples"] plot_chains(radiosamples, errors[mask_radio_pars],labels[mask_radio_pars], savepath="{0:s}_{1:s}.png".format(options.resultsfile,"samples_radio"),debug=options.debug) if dogamma+doboth == 1 or options.compare == 1: gammalabels = labels[mask_gamma_pars] if options.plotonly==0: if doboth == 0: startpos = np.copy(pos) else: startpos = np.copy(jragpos[:,mask_gamma_pars]) gammasampler = emcee.EnsembleSampler( nwalker,len(mask_gamma_pars),log_prob_quick, #args=(labels, values, errors, npeaks,M_gamma, phases,weights,M_radio,radio_residuals,radio_toa_uncerts), args =(gammalabels,values[mask_gamma_pars],errors[mask_gamma_pars],npeaks,M_gamma[:,mask_gamma_pars[:-3*npeaks]],phases,weights,None, None, None), moves=[(emcee.moves.KDEMove(), 1.0),],) if nburn > 0: gammasampler,gburnpos = mcmc_run(gammasampler,startpos,nburn,usetau=0,doburnin=1) gammasamples = gammasampler.get_chain() plot_chains(gammasamples, errors[mask_gamma_pars],labels[mask_gamma_pars], savepath="{0:s}_{1:s}.png".format(options.resultsfile,"burnin_gamma"),debug=options.debug) else: gburnpos = np.copy(startpos) gammasampler.reset() gammasampler,gammapos,ngthin = mcmc_run(gammasampler,gburnpos,nsteps,usetau=options.usetau,doburnin=0, savepath="{0:s}_{1:s}.png".format(options.resultsfile,"converg_gamma"),debug=options.debug) gammasamples = gammasampler.get_chain() if options.usetau: flat_gammasamples = gammasampler.get_chain(discard=20*ngthin, thin=ngthin, flat=True) log_prob_gammasamples = gammasampler.get_log_prob(discard=20*ngthin, thin=ngthin, flat=True) else: flat_gammasamples = gammasampler.get_chain(discard=0, thin=ngthin, flat=True) log_prob_gammasamples = gammasampler.get_log_prob(discard=0, thin=ngthin, flat=True) # store samples in compressed file gammasamplesdict = {} gammasamplesdict["gammasamples"] = gammasamples gammasamplesdict["flat_gammasamples"] = flat_gammasamples gammasamplesdict["log_prob_gammasamples"] = log_prob_gammasamples np.savez_compressed("{0:s}_chain_gamma".format(options.resultsfile), **gammasamplesdict) loaded = np.load("{0:s}_chain_gamma.npz".format(options.resultsfile)) gammasamples = loaded["gammasamples"] flat_gammasamples = loaded["flat_gammasamples"] log_prob_gammasamples = loaded["log_prob_gammasamples"] plot_chains(gammasamples, errors[mask_gamma_pars],labels[mask_gamma_pars], savepath="{0:s}_{1:s}.png".format(options.resultsfile,"samples_gamma"),debug=options.debug) ############################################################################ # MAIN: Visualization of MCMC results ############################################################################ if doboth+doradio == 1: flat_samples = np.copy(flat_radiosamples) log_prob_samples = np.copy(log_prob_radiosamples) elif doboth+dogamma == 1: flat_samples = np.copy(flat_gammasamples) log_prob_samples = np.copy(log_prob_gammasamples) maxL_theta = flat_samples[np.where(log_prob_samples == np.max(log_prob_samples))][0] maxL_val = values+errors*maxL_theta if options.debug > 2: print("\n# Values, maxL-values, shift") print(values) print(maxL_val) print(np.array(values)-np.array(maxL_val)) cen_val = np.zeros(nfitpars) upp_val = np.zeros(nfitpars) low_val = np.zeros(nfitpars) err_val = np.zeros(nfitpars) err_val_max = np.zeros(nfitpars) for i in range(nfitpars): low_val[i],cen_val[i],upp_val[i] = values[i]+errors[i]*np.percentile(flat_samples[:, i],[16,50,84]) err_val_max[i] = np.max([maxL_val[i]-low_val[i],upp_val[i]-maxL_val[i]]) err_val[i] = errors[i]*np.std(flat_samples[:,i]) if options.debug > 1: print("{0} {1} {2} {3}".format(low_val[i],cen_val[i],upp_val[i],maxL_val[i])) print("{0} {1} {2}".format(errors[i]*np.std(flat_samples[:, i]),maxL_val[i]-low_val[i],upp_val[i]-maxL_val[i])) print(err_val_max[i]) print(err_val[i]) print("#####################") if "EPS1" in labels: print("95 percent upper limit on ECC") id_par1 = np.where(labels == "EPS1")[0][0] par1_samples = values[id_par1]+errors[id_par1]*flat_samples[:, id_par1] print("EPS1 ",np.percentile(par1_samples,[5,95]),maxL_val[id_par1],np.percentile(par1_samples,[5,95])-maxL_val[id_par1]) id_par2 = np.where(labels == "EPS2")[0][0] par2_samples = values[id_par2]+errors[id_par2]*flat_samples[:, id_par2] print("EPS2 ",np.percentile(par2_samples,[5,95]),maxL_val[id_par2],np.percentile(par2_samples,[5,95])-maxL_val[id_par2]) print("ECC ",np.percentile(np.sqrt(par1_samples**2+par2_samples**2),[5,95]), np.sqrt(maxL_val[id_par1]**2+maxL_val[id_par2]**2), "-",np.sqrt(maxL_val[id_par1]**2+maxL_val[id_par2]**2)-np.percentile(np.sqrt(par1_samples**2+par2_samples**2),[16]), "+",np.percentile(np.sqrt(par1_samples**2+par2_samples**2),[84])-np.sqrt(maxL_val[id_par1]**2+maxL_val[id_par2]**2)) print("OM ",np.arctan(maxL_val[id_par1]/maxL_val[id_par2])/np.pi*180.0, (np.arctan(maxL_val[id_par1]/maxL_val[id_par2])-np.percentile(np.arctan(par1_samples/par2_samples),[16,84]))/np.pi*180.0) if "PMRA" in labels: print("95 percent upper limit on PM") id_par1 = np.where(labels == "PMRA")[0][0] par1_samples = values[id_par1]+errors[id_par1]*flat_samples[:, id_par1] print("PMRA ",np.percentile(par1_samples,[5,95]),maxL_val[id_par1],np.percentile(par1_samples,[5,95])-maxL_val[id_par1]) id_par2 = np.where(labels == "PMDEC")[0][0] par2_samples = values[id_par2]+errors[id_par2]*flat_samples[:, id_par2] print("PMDEC ",np.percentile(par2_samples,[5,95]),maxL_val[id_par2],np.percentile(par2_samples,[5,95])-maxL_val[id_par2]) print("PM ",np.percentile(np.sqrt(par1_samples**2+par2_samples**2),[5,95]), np.sqrt(maxL_val[id_par1]**2+maxL_val[id_par2]**2), "-",np.sqrt(maxL_val[id_par1]**2+maxL_val[id_par2]**2)-np.percentile(np.sqrt(par1_samples**2+par2_samples**2),[16]), "+",np.percentile(np.sqrt(par1_samples**2+par2_samples**2),[84])-np.sqrt(maxL_val[id_par1]**2+maxL_val[id_par2]**2)) if "F2" in labels: print("95 percent upper limit on F2") id_par1 = np.where(labels == "F2")[0][0] par1_samples = values[id_par1]+errors[id_par1]*flat_samples[:, id_par1] print("F2 ",np.percentile(par1_samples,[5,95]),maxL_val[id_par1],np.percentile(par1_samples,[5,95])-maxL_val[id_par1]) if "FB1" in labels: print("95 percent upper limit on FB1") id_par1 = np.where(labels == "FB1")[0][0] par1_samples = values[id_par1]+errors[id_par1]*flat_samples[:, id_par1] print("FB1 ",np.percentile(par1_samples,[5,95]),maxL_val[id_par1],np.percentile(par1_samples,[5,95])-maxL_val[id_par1]) print("#####################") nbins = 50 CORNER_KWARGS = dict( plot_datapoints=False, plot_density=True, fill_contours=True, plot_contours=True, levels=(0.0, 1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.)), bins=nbins, ) if doboth: bdeltas = (values+errors*flat_samples) - maxL_val bdelta_maxs = np.max(np.abs(bdeltas),axis=0) bdelta_mags = order_of_mag(bdelta_maxs) bdelta_pres = bdeltas/(10.0**bdelta_mags) new_blabels = [] for i in range(len(labels)): if junits[i] == "1": new_blabels.append(r"$\Delta$ {0:s} [$10^{1:s}{2:d}{3:s}$]".format(labels[i],"{",bdelta_mags[i],"}")) else: new_blabels.append(r"$\Delta$ {0:s} [$10^{1:s}{2:d}{3:s}$ {4:s}]".format(labels[i],"{",bdelta_mags[i],"}",junits[i])) figure = corner.corner(bdelta_pres, labels=new_blabels,truths=np.zeros_like(maxL_val),**CORNER_KWARGS, ); if options.resultsfile is not None: plt.savefig("{0:s}_triplot.png".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() bcovmat = np.cov(flat_samples,rowvar=False) if options.debug > 0: print("Covariance matrix for joint timing analysis:") print(bcovmat) print(np.shape(bcovmat)) np.savetxt("{0:s}_covmat.cov".format(options.resultsfile),bcovmat) if doradio+doboth == 1 or options.compare == 1:#doradio: rdeltas = (values[mask_radio_pars]+errors[mask_radio_pars]*flat_radiosamples) - maxL_val[mask_radio_pars] rdelta_maxs = np.max(np.abs(rdeltas),axis=0) rdelta_mags = order_of_mag(rdelta_maxs) if doboth == 0: rdelta_pres = rdeltas/(10.0**rdelta_mags) new_rlabels = [] for i in range(len(labels)): if junits[i] == "1": new_rlabels.append(r"$\Delta$ {0:s} [$10^{1:s}{2:d}{3:s}$]".format(labels[i],"{",rdelta_mags[i],"}")) else: new_rlabels.append(r"$\Delta$ {0:s} [$10^{1:s}{2:d}{3:s}$ {4:s}]".format(labels[i],"{",rdelta_mags[i],"}",junits[i])) figure = corner.corner(rdelta_pres, labels=new_rlabels,truths=np.zeros_like(maxL_val),color="r",**CORNER_KWARGS, ); if options.resultsfile is not None: plt.savefig("{0:s}_triplot_radio.png".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() elif options.compare: figure = corner.corner(rdeltas/(10.0**bdelta_mags[mask_radio_pars]), **CORNER_KWARGS,labels=np.array(new_blabels)[mask_radio_pars],color="r", hist_kwargs=dict(histtype='stepfilled',ec=mpl.colors.to_rgba("r",alpha=None),fc=mpl.colors.to_rgba("r",alpha=0.4),lw=1.5), weights=np.ones(len(flat_radiosamples))*len(flat_samples[:,mask_radio_pars])/len(flat_radiosamples) ); corner.corner(bdelta_pres[:,mask_radio_pars],fig=figure, **CORNER_KWARGS,labels=np.array(new_blabels)[mask_radio_pars],color="k", hist_kwargs=dict(histtype='stepfilled',ec=mpl.colors.to_rgba("k",alpha=None),fc=mpl.colors.to_rgba("k",alpha=0.4),lw=1.5), ); if options.resultsfile is not None: plt.savefig("{0:s}_triplot_radio.png".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() if dogamma+doboth == 1 or options.compare == 1:#dogamma: gdeltas = (values[mask_gamma_pars]+errors[mask_gamma_pars]*flat_gammasamples) - maxL_val[mask_gamma_pars] gdelta_maxs = np.max(np.abs(gdeltas),axis=0) gdelta_mags = order_of_mag(gdelta_maxs) if doboth == 0: gdelta_pres = gdeltas/(10.0**gdelta_mags) new_glabels = [] for i in range(len(labels)): if junits[i] == "1": new_glabels.append(r"$\Delta$ {0:s} [$10^{1:s}{2:d}{3:s}$]".format(labels[i],"{",gdelta_mags[i],"}")) else: new_glabels.append(r"$\Delta$ {0:s} [$10^{1:s}{2:d}{3:s}$ {4:s}]".format(labels[i],"{",gdelta_mags[i],"}",junits[i])) figure = corner.corner(gdelta_pres, labels=new_glabels,truths=np.zeros_like(maxL_val),color="b",**CORNER_KWARGS, ); if options.resultsfile is not None: plt.savefig("{0:s}_triplot_gamma.png".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() elif options.compare: figure = corner.corner(gdeltas/(10.0**bdelta_mags[mask_gamma_pars]), **CORNER_KWARGS,labels=np.array(new_blabels)[mask_gamma_pars],color="b", hist_kwargs=dict(histtype='stepfilled',ec=mpl.colors.to_rgba("b",alpha=None),fc=mpl.colors.to_rgba("b",alpha=0.4),lw=1.5), weights=np.ones(len(flat_gammasamples))*len(flat_samples[:,mask_gamma_pars])/len(flat_gammasamples), ); corner.corner(bdelta_pres[:,mask_gamma_pars],fig=figure, **CORNER_KWARGS,labels=np.array(new_blabels)[mask_gamma_pars],color="k", hist_kwargs=dict(histtype='stepfilled',ec=mpl.colors.to_rgba("k",alpha=None),fc=mpl.colors.to_rgba("k",alpha=0.4),lw=1.5), ); if options.resultsfile is not None: plt.savefig("{0:s}_triplot_gamma.png".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() # -------------------------------------------------------------------------- if options.compare == 1: fig, ax = plt.subplots(ncols=nfitpars,nrows=nfitpars,figsize=(15,15)) cornerrange=[] for n in range(nfitpars): cornerrange.append((np.min(np.concatenate(((flat_samples[:,np.where(labels==labels[n])]).flatten(), (flat_radiosamples[:,np.where(radiolabels==labels[n])]).flatten(), (flat_gammasamples[:,np.where(gammalabels==labels[n])]).flatten()))), np.max(np.concatenate(((flat_samples[:,np.where(labels==labels[n])]).flatten(), (flat_radiosamples[:,np.where(radiolabels==labels[n])]).flatten(), (flat_gammasamples[:,np.where(gammalabels==labels[n])]).flatten()))))) cornerrange = np.array(cornerrange) for r in range(nfitpars): for c in range(nfitpars): if c > r: ax[r,c].axis("off") elif c == r: ax[c,c].set_yticklabels([]) if labels[c] in radiolabels: c_tmp = np.where(radiolabels == labels[c]) ax[c,c].hist(values[c]+errors[c]*(flat_radiosamples[:,c_tmp]).flatten(), bins=nbins,#range=(values[c]+errors[c]*cornerrange[c]), histtype="stepfilled",density=True, ec=mpl.colors.to_rgba("r",alpha=None),fc=mpl.colors.to_rgba("r",alpha=0.4),lw=1.5, ); if labels[c] in gammalabels: c_tmp = np.where(gammalabels == labels[c]) ax[c,c].hist(values[c]+errors[c]*(flat_gammasamples[:,c_tmp]).flatten(), bins=nbins,#range=(values[c]+errors[c]*cornerrange[c]), histtype="stepfilled",density=True, ec=mpl.colors.to_rgba("b",alpha=None),fc=mpl.colors.to_rgba("b",alpha=0.4),lw=1.5, ); ax[c,c].hist(values[c]+errors[c]*(flat_samples[:,c]).flatten(), bins=nbins,#range=(values[c]+errors[c]*cornerrange[c]), histtype="stepfilled",density=True, ec=mpl.colors.to_rgba("k",alpha=None),fc=mpl.colors.to_rgba("k",alpha=0.4),lw=1.5, ); if r == nfitpars-1: ax[c,c].set_xlabel(labels[c]) ax[c,c].set_xticklabels(ax[c,c].get_xticklabels(),rotation=45) else: ax[c,c].set_xticklabels([]) else: if labels[c] in radiolabels and labels[r] in radiolabels: c_tmp = np.where(radiolabels == labels[c]) r_tmp = np.where(radiolabels == labels[r]) corner.hist2d(values[c]+errors[c]*flat_radiosamples[:,c_tmp], values[r]+errors[r]*flat_radiosamples[:,r_tmp], ax=ax[r,c],**CORNER_KWARGS,color="r", range=(values[c]+errors[c]*cornerrange[c],values[r]+errors[r]*cornerrange[r]), ); if labels[c] in gammalabels and labels[r] in gammalabels: c_tmp = np.where(gammalabels == labels[c]) r_tmp = np.where(gammalabels == labels[r]) corner.hist2d(values[c]+errors[c]*flat_gammasamples[:,c_tmp], values[r]+errors[r]*flat_gammasamples[:,r_tmp], ax=ax[r,c],**CORNER_KWARGS,color="b", range=(values[c]+errors[c]*cornerrange[c],values[r]+errors[r]*cornerrange[r]), ); corner.hist2d(values[c]+errors[c]*flat_samples[:,c], values[r]+errors[r]*flat_samples[:,r], ax=ax[r,c],**CORNER_KWARGS,color="k", range=(values[c]+errors[c]*cornerrange[c],values[r]+errors[r]*cornerrange[r]) ); if c == 0: ax[r,c].set_ylabel(labels[r]) ax[r,c].set_yticklabels(ax[r,c].get_yticklabels(),rotation=45) else: ax[r,c].set_yticklabels([]) if r == nfitpars-1: ax[r,c].set_xlabel(labels[c]) ax[r,c].set_xticklabels(ax[r,c].get_xticklabels(),rotation=45) else: ax[r,c].set_xticklabels([]) plt.subplots_adjust(wspace=0.03, hspace=0.03) if options.resultsfile is not None: plt.savefig("{0:s}_triplot_compare_all.png".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() fig, ax = plt.subplots(ncols=len(compareskylabels),nrows=len(compareskylabels),figsize=(15,15)) cornerrange=[] for n in range(len(compareskylabels)): cornerrange.append((np.min(np.concatenate(((flat_samples[:,np.where(labels==compareskylabels[n])]).flatten(), (flat_radiosamples[:,np.where(radiolabels==compareskylabels[n])]).flatten(), (flat_gammasamples[:,np.where(gammalabels==compareskylabels[n])]).flatten()))), np.max(np.concatenate(((flat_samples[:,np.where(labels==compareskylabels[n])]).flatten(), (flat_radiosamples[:,np.where(radiolabels==compareskylabels[n])]).flatten(), (flat_gammasamples[:,np.where(gammalabels==compareskylabels[n])]).flatten()))))) cornerrange = np.array(cornerrange) for r in range(len(compareskylabels)): for c in range(len(compareskylabels)): if c > r: ax[r,c].axis("off") elif c == r: ax[c,c].set_yticklabels([]) c_main = np.where(labels == compareskylabels[c]) if compareskylabels[c] in radiolabels: c_tmp = np.where(radiolabels == compareskylabels[c]) ax[c,c].hist(values[c_main]+errors[c_main]*(flat_radiosamples[:,c_tmp]).flatten(), bins=nbins,range=(values[c_main]+errors[c_main]*cornerrange[c]), histtype="stepfilled",density=True, ec=mpl.colors.to_rgba("r",alpha=None),fc=mpl.colors.to_rgba("r",alpha=0.4),lw=1.5, ); if compareskylabels[c] in gammalabels: c_tmp = np.where(gammalabels == compareskylabels[c]) ax[c,c].hist(values[c_main]+errors[c_main]*(flat_gammasamples[:,c_tmp]).flatten(), bins=nbins,range=(values[c_main]+errors[c_main]*cornerrange[c]), histtype="stepfilled",density=True, ec=mpl.colors.to_rgba("b",alpha=None),fc=mpl.colors.to_rgba("b",alpha=0.4),lw=1.5, ); ax[c,c].hist(values[c_main]+errors[c_main]*(flat_samples[:,c_tmp]).flatten(), bins=nbins,range=(values[c_main]+errors[c_main]*cornerrange[c]), histtype="stepfilled",density=True, ec=mpl.colors.to_rgba("k",alpha=None),fc=mpl.colors.to_rgba("k",alpha=0.4),lw=1.5, ); if r == len(compareskylabels)-1: ax[c,c].set_xlabel(compareskylabels[c]) ax[c,c].set_xticklabels(ax[c,c].get_xticklabels(),rotation=45) else: ax[c,c].set_xticklabels([]) else: c_main = np.where(labels == compareskylabels[c]) r_main = np.where(labels == compareskylabels[r]) if compareskylabels[c] in radiolabels and compareskylabels[r] in radiolabels: c_tmp = np.where(radiolabels == compareskylabels[c]) r_tmp = np.where(radiolabels == compareskylabels[r]) corner.hist2d(values[c_main]+errors[c_main]*flat_radiosamples[:,c_tmp], values[r_main]+errors[r_main]*flat_radiosamples[:,r_tmp], ax=ax[r,c],**CORNER_KWARGS,color="r", range=(values[c_main]+errors[c_main]*cornerrange[c],values[r_main]+errors[r_main]*cornerrange[r]), ); if compareskylabels[c] in gammalabels and compareskylabels[r] in gammalabels: c_tmp = np.where(gammalabels == compareskylabels[c]) r_tmp = np.where(gammalabels == compareskylabels[r]) corner.hist2d(values[c_main]+errors[c_main]*flat_gammasamples[:,c_tmp], values[r_main]+errors[r_main]*flat_gammasamples[:,r_tmp], ax=ax[r,c],**CORNER_KWARGS,color="b", range=(values[c_main]+errors[c_main]*cornerrange[c],values[r_main]+errors[r_main]*cornerrange[r]), ); corner.hist2d(values[c_main]+errors[c_main]*flat_samples[:,c_main], values[r_main]+errors[r_main]*flat_samples[:,r_main], ax=ax[r,c],**CORNER_KWARGS,color="k", range=(values[c_main]+errors[c_main]*cornerrange[c],values[r_main]+errors[r_main]*cornerrange[r]), ); if c == 0: ax[r,c].set_ylabel(compareskylabels[r]) ax[r,c].set_yticklabels(ax[r,c].get_yticklabels(),rotation=45) else: ax[r,c].set_yticklabels([]) if r == len(compareskylabels)-1: ax[r,c].set_xlabel(compareskylabels[c]) ax[r,c].set_xticklabels(ax[r,c].get_xticklabels(),rotation=45) else: ax[r,c].set_xticklabels([]) plt.subplots_adjust(wspace=0.03, hspace=0.03) if options.resultsfile is not None: plt.savefig("{0:s}_triplot_compare.png".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() # -------- figure = corner.corner(rdeltas[:,mask_compare_radio]/(10.0**bdelta_mags[mask_compare_both]), **CORNER_KWARGS,labels=np.array(new_blabels)[mask_compare_both],color="r", hist_kwargs=dict(histtype='stepfilled',ec=mpl.colors.to_rgba("r",alpha=None),fc=mpl.colors.to_rgba("r",alpha=0.4),lw=1.5), ); corner.corner(gdeltas[:,mask_compare_gamma]/(10.0**bdelta_mags[mask_compare_both]),fig=figure, **CORNER_KWARGS,labels=np.array(new_blabels)[mask_compare_both],color="b", hist_kwargs=dict(histtype='stepfilled',ec=mpl.colors.to_rgba("b",alpha=None),fc=mpl.colors.to_rgba("b",alpha=0.4),lw=1.5), ); corner.corner(bdelta_pres[:,mask_compare_both],fig=figure, **CORNER_KWARGS,labels=np.array(new_blabels)[mask_compare_both],color="k", hist_kwargs=dict(histtype='stepfilled',ec=mpl.colors.to_rgba("k",alpha=None),fc=mpl.colors.to_rgba("k",alpha=0.4),lw=1.5), ); if options.resultsfile is not None: plt.savefig("{0:s}_triplot_compare2.pdf".format(options.resultsfile),bbox_inches='tight') if options.debug > 1: plt.show() else: plt.show() ################################################################################ # MAIN: Report results ################################################################################ print(fitpars) for i in range(nfitpars): if labels[i] == "F1": print("{0:s} {1:.14g}".format(labels[i],maxL_val[i])) else: print("{0:s} {1:.14f}".format(labels[i],maxL_val[i])) print("\nShift as in MCMC") if doboth: print("Using fast function: {0: 14.2f}".format(log_prob_quick(maxL_theta,labels,values,errors,npeaks,M_gamma,phases,weights,M_radio,radio_residuals,radio_toa_uncerts))) if dogamma: print("gamma only: {0: 14.2f}".format(log_prob_quick(maxL_theta,labels,values,errors,npeaks,M_gamma,phases,weights,None,None,None))) if doradio: print("radio only: {0: 14.2f}".format(log_prob_quick(maxL_theta,labels,values,errors,npeaks,None,None,None,M_radio,radio_residuals,radio_toa_uncerts))) if doboth == 1: for i,p in enumerate(fitpars): getattr(timing_model,p).value = maxL_val[:-npeaks*3][i] elif dogamma == 1: if len(jumppars)>0: for i,p in enumerate(fitpars[:-len(jumppars)]): getattr(timing_model,p).value = maxL_val[:-npeaks*3][i] else: for i,p in enumerate(fitpars): getattr(timing_model,p).value = maxL_val[:-npeaks*3][i] elif doradio == 1: for i,p in enumerate(fitpars): getattr(timing_model,p).value = maxL_val[i] if dogamma: phases = np.mod(timing_model.phase(fermi_toas).frac,1.0) latprofile_new = templ_gen(maxL_val[-npeaks*3:],npeaks) if doradio: radio_residuals = Residuals(radio_toas,timing_model).phase_resids maxL_realval = np.copy(maxL_val) theta_null = np.zeros(nfitpars) print("\nFolding with new pars") print(fitpars) if doboth: print("Using function: + {0: 14.2f}".format(log_prob_precise(theta_null,fitpars,labels,maxL_val,errors,npeaks,fermi_toas,weights,radio_toas,radio_toa_uncerts,timing_model))) if dogamma: if len(jumppars)>0: print("gamma only: {0: 14.2f}".format(log_prob_precise(theta_null,fitpars[:-len(fitpars)],labels,maxL_val,errors,npeaks,fermi_toas,weights,None,None,timing_model))) else: print("gamma only: {0: 14.2f}".format(log_prob_precise(theta_null,fitpars,labels,maxL_val,errors,npeaks,fermi_toas,weights,None,None,timing_model))) if doradio: print("radio only: {0: 14.2f}".format(log_prob_precise(theta_null,fitpars,labels,maxL_val,errors,npeaks,None,None,radio_toas,radio_toa_uncerts,timing_model))) # update of Mgamma and Mradio not required due to theta = 0 if doboth: print("Using fast function: {0: 14.2f}".format(log_prob_quick(theta_null,labels,maxL_val,errors,npeaks,M_gamma,phases,weights,M_radio,radio_residuals,radio_toa_uncerts))) if dogamma: print("gamma only: {0: 14.2f}".format(log_prob_quick(theta_null,labels,maxL_val,errors,npeaks,M_gamma,phases,weights,None,None,None))) if doradio: print("radio only: {0: 14.2f}".format(log_prob_quick(theta_null,labels,maxL_val,errors,npeaks,None,None,None,M_radio,radio_residuals,radio_toa_uncerts))) if doboth: plot_phases(fermi_toas,weights,radio_toas,timing_model,template=latprofile_new,isolated=isolated,savepath="{0:s}_phases.pdf".format(options.resultsfile)) elif doradio: plot_radio_residuals(radio_toas,timing_model,isolated=isolated,savepath="{0:s}_radiopostfit.png".format(options.resultsfile),debug=options.debug) ################################################################################ # MAIN: Store results ################################################################################ # store LAT photons in tvsph file if dogamma: fermi_t = timing_model.get_barycentric_toas(fermi_toas).value start = np.floor(np.min(fermi_t)) finish = np.ceil(np.max(fermi_t)) fermi_orbphi = timing_model.orbital_phase(fermi_toas,anom="mean",radians=False) with open("{0:s}.tvsph".format(options.resultsfile),"w") as tvsph_file: for i in range(len(fermi_t)): if isolated == False: tvsph_file.write("{0:.17f} {1:.15f} {2:.6f} {3:.5f} {4:.6f}\n".format(fermi_t[i],phases[i],weights[i],energies[i],fermi_orbphi[i])) else: tvsph_file.write("{0:.17f} {1:.15f} {2:.6f} {3:.5f}\n".format(fermi_t[i],phases[i],weights[i],energies[i])) elif doradio: radio_t = timing_model.get_barycentric_toas(radio_toas).value start = np.floor(np.min(radio_t)) finish = np.ceil(np.max(radio_t)) if dogamma: # store best-fit tempfile with open("{0:s}.itemp".format(options.resultsfile),"w") as ifile: ifile.write("# gauss\n") ifile.write("-------------------------\n") const_value = 1.0 const_error = 0.0 for i in range(npeaks): if i == 0: id_par = np.where(labels == "PHI0")[0][0] tmp_phi0 = maxL_realval[id_par] parvalue = np.copy(tmp_phi0) else: id_par = np.where(labels == "DPHI{0:d}".format(i))[0][0] parvalue = np.mod(maxL_realval[id_par]+tmp_phi0,1.0) parerror = np.copy(err_val[id_par]) ifile.write("phas{0:d} = {1:.8f} +/- {2:.8f}\n".format(i+1,parvalue,parerror)) id_par = np.where(labels == "FWHM{0:d}".format(i))[0][0] parvalue = np.copy(maxL_realval[id_par]) parerror = np.copy(err_val[id_par]) ifile.write("fwhm{0:d} = {1:.8f} +/- {2:.8f}\n".format(i+1,parvalue,parerror)) id_par = np.where(labels == "AMPL{0:d}".format(i))[0][0] parvalue = np.copy(maxL_realval[id_par]) const_value -= parvalue parerror = np.copy(err_val[id_par]) const_error += parerror ifile.write("ampl{0:d} = {1:.8f} +/- {2:.8f}\n".format(i+1,parvalue,parerror)) ifile.write("const = {0:.8f} +/- {1:.8f}\n".format(const_value,const_error)) ifile.write("-------------------------\n") # store best-fit parfile with open(options.parfile,"r") as pfile: lines = pfile.readlines() with open("{0:s}.par".format(options.resultsfile),"w") as nfile: writtenstart,writtenfinish = 0,0 print("\n### New parfile : ###\n") for line in lines: splittedline = line.strip().split() linelen = len(splittedline) if linelen == 0: continue if line[0] == "#": continue lineparts = np.copy(splittedline).tolist() if splittedline[0] in fitpars: if splittedline[2] == "1": id_par = np.where(labels == splittedline[0])[0][0] if splittedline[0] == "RAJ": raj = Angle(maxL_realval[id_par], u.hour) parvalue = "{0:02.0f}:".format(np.abs(raj.hms.h)) parvalue += "{0:02.0f}:".format(np.abs(raj.hms.m)) if np.abs(raj.hms.s) < 10.0: parvalue += "0{0:.15f}".format(np.abs(raj.hms.s)) else: parvalue += "{0:.15f}".format(np.abs(raj.hms.s)) parerror = "{0:.14g}".format(err_val[id_par]*3600.0) elif splittedline[0] == "DECJ": decj = Angle(maxL_realval[id_par], u.degree) if decj.dms.d < 0.0: parvalue = "-" else: parvalue = "+" parvalue += "{0:02.0f}:".format(np.abs(decj.dms.d)) parvalue += "{0:02.0f}:".format(np.abs(decj.dms.m)) if np.abs(decj.dms.s) < 10.0: parvalue += "0{0:.15f}".format(np.abs(decj.dms.s)) else: parvalue += "{0:.15f}".format(np.abs(decj.dms.s)) parerror = "{0:.14g}".format(err_val[id_par]*3600.0) else: parvalue = "{0:.16g}".format(maxL_realval[id_par]) parerror = "{0:.14g}".format(err_val[id_par]) lineparts[1] = parvalue if linelen > 4: lineparts[4] = parerror else: lineparts.append(parerror) linelen += 1 if splittedline[0] == "START": lineparts[1] = start writtenstart = 1 if splittedline[0] == "FINISH": lineparts[1] = finish writtenfinish = 1 if "JUMP" in splittedline[0]: if splittedline[2] in jumptels and linelen == 6: # uncertainty given, but could be 0 if doradio: # else: no search in any jump par id_jump = np.where(np.array(jumptels) == splittedline[2])[0][0] id_par = np.where(labels == jumppars[id_jump])[0][0] parvalue = "{0:.16g}".format(maxL_realval[id_par]) lineparts[3] = parvalue if lineparts[4] == "1": parerror = "{0:.14g}".format(err_val[id_par]) lineparts[5] = parerror if lineparts[3][0] in ["+","-"]: newline = "JUMP {0} {1:<10} {2:<28} {3:<4} {4}\n".format(lineparts[1],lineparts[2],lineparts[3],lineparts[4],lineparts[5]) else: newline = "JUMP {0} {1:<11} {2:<27} {3:<4} {4}\n".format(lineparts[1],lineparts[2],lineparts[3],lineparts[4],lineparts[5]) elif linelen == 6: # no search in this jump par if lineparts[3][0] in ["+","-"]: newline = "JUMP {0} {1:<10} {2:<28} {3:<4} {4}\n".format(lineparts[1],lineparts[2],lineparts[3],lineparts[4],lineparts[5]) else: newline = "JUMP {0} {1:<11} {2:<27} {3:<4} {4}\n".format(lineparts[1],lineparts[2],lineparts[3],lineparts[4],lineparts[5]) elif linelen == 5: # no search in this jump par if lineparts[3][0] in ["+","-"]: newline = "JUMP {0} {1:<10} {2:<28} {3}\n".format(lineparts[1],lineparts[2],lineparts[3],lineparts[4]) else: newline = "JUMP {0} {1:<11} {2:<27} {3}\n".format(lineparts[1],lineparts[2],lineparts[3],lineparts[4]) elif linelen > 2: if lineparts[1][0] in ["+","-"]: newline = "{0:<18} {1:<29}".format(lineparts[0],lineparts[1]) else: newline = "{0:<19} {1:<28}".format(lineparts[0],lineparts[1]) if linelen == 3: newline += "{0}\n".format(lineparts[2]) elif linelen > 3: newline += "{0:<5}".format(lineparts[2]) if linelen == 4: newline += "{0}\n".format(lineparts[3]) elif linelen > 4: newline += "{0:<10}".format(lineparts[3]) newline += "{0}\n".format(lineparts[4]) elif linelen == 2: newline = "{0:<19} {1}\n".format(lineparts[0],lineparts[1]) nfile.write(newline) print(newline.strip()) if writtenstart == 0: newline = "{0:<19} {1}\n".format("START",start) nfile.write(newline) print(newline.strip()) if writtenfinish == 0: newline = "{0:<19} {1}\n".format("FINISH",finish) nfile.write(newline) print(newline.strip()) ################################################################################ # MAIN: Prepare for next run ################################################################################ if dogamma: shutil.copy("{0:s}.itemp".format(options.resultsfile),options.newtempfile) shutil.copy("{0:s}.par".format(options.resultsfile),options.newparfile)