From c3693707d9c91c01c4a523edfb057c251c90fe40 Mon Sep 17 00:00:00 2001 From: wbalmer Date: Tue, 2 Jan 2024 14:07:31 -0500 Subject: [PATCH] create reference epoch, stationary case, and chi2 objects. refigured ref_epoch --- .gitignore | 6 +- backtracks/backtracks.py | 88 ++++++++++++++++---- backtracks/plotting.py | 171 ++++++++++++++++++++++++++++----------- 3 files changed, 203 insertions(+), 62 deletions(-) diff --git a/.gitignore b/.gitignore index cfeb488..3336246 100644 --- a/.gitignore +++ b/.gitignore @@ -14,10 +14,12 @@ bgstar_emcee.py .vscode* *beta_pic* *SAO* +*sao* +*yses2* +*YSES2* .pyproject.toml.sw .pytest_cache/ +tests/test_changes.py tests/HD_131399A_* tests/dynesty.save -tests/yses2* -tests/YSES_2* tests/gaia_query_* \ No newline at end of file diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index 11da571..73b75d2 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -22,7 +22,7 @@ import numpy as np import pandas as pd -from astropy.coordinates import SkyCoord +from astropy.coordinates import SkyCoord, Distance from astropy.table import Table from astropy.io import fits from astropy.time import Time @@ -34,7 +34,7 @@ from schwimmbad import MPIPool from backtracks.utils import pol2car, transform_gengamm, transform_normal, transform_uniform, utc2tt, HostStarPriors -from backtracks.plotting import diagnostic, neighborhood, plx_prior, posterior, trackplot +from backtracks.plotting import diagnostic, neighborhood, plx_prior, posterior, trackplot, stationtrackplot # Set the Gaia data release to DR3 @@ -162,6 +162,7 @@ def __init__(self, target_name: str, candidate_file: str, nearby_window: float = if kwargs['rv_host_method'].lower() == 'uniform': self.rv_host_method='uniform' self.rv_lower,self.rv_upper=kwargs['rv_host_params'] # e.g., rv_host_params=(-10,10) + self.radvelo = np.mean([self.rv_lower,self.rv_upper]) elif kwargs['rv_host_method'].lower() == 'normal': self.rv_host_method='normal' self.radvelo,self.sig_rv=kwargs['rv_host_params'] # e.g., rv_host_params=(10,0.5) @@ -180,17 +181,40 @@ def __init__(self, target_name: str, candidate_file: str, nearby_window: float = [ra_dec_corr*sig_ra*sig_dec,sig_dec**2,sig_dec*sig_pmra*dec_pmra_corr,sig_dec*sig_pmdec*dec_pmdec_corr,sig_dec*sig_parallax*dec_parallax_corr], [ra_pmra_corr*sig_ra*sig_pmra,dec_pmra_corr*sig_pmra*sig_dec,sig_pmra**2,sig_pmra*sig_pmdec*pmra_pmdec_corr,sig_pmra*sig_parallax*parallax_pmra_corr],[sig_pmdec*sig_ra*ra_pmdec_corr,sig_pmdec*sig_dec*dec_pmdec_corr,sig_pmdec*sig_pmra*pmra_pmdec_corr,sig_pmdec**2,sig_pmdec*sig_parallax*parallax_pmdec_corr],[sig_parallax*sig_ra*ra_parallax_corr,sig_parallax*sig_dec*dec_parallax_corr,sig_parallax*sig_pmra*parallax_pmra_corr,sig_parallax*sig_pmdec*parallax_pmdec_corr,sig_parallax**2]]) self.HostStarPriors = HostStarPriors(self.host_mean,self.host_cov) + # initial estimate for background star scenario (best guesses) - # (manually looked at approximate offset from star with cosine compensation for - # dependence of RA on declination (RA is a singularity at the declination poles)) - # this offset is calculated wrong (should use spherical offsets to) but it kinda works because the uniform prior is bigger than this. not all of these parameters are still in use - self.ra0 = self.rao-(self.ras[0])/1000/3600/np.cos(self.deco/180.*np.pi) - self.dec0 = self.deco-(self.decs[0])/1000/3600 + init_host_coord = SkyCoord(ra=self.rao*u.deg, dec=self.deco*u.deg, distance=Distance(parallax=self.paro*u.mas), pm_ra_cosdec=self.pmrao*u.mas/u.yr, pm_dec=self.pmdeco*u.mas/u.yr, obstime=Time(self.gaia_epoch,format='jyear',scale='tcb')) + # choose time of reference observation of relative candidate position (defaults to the first observation) + if 'ref_epoch_idx' in kwargs: + self.ref_epoch_idx = kwargs['ref_epoch_idx'] + self.ref_epoch = self.epochs[self.ref_epoch_idx] + else: + self.ref_epoch_idx = 0 + self.ref_epoch = self.epochs[self.ref_epoch_idx] + # propagate host position from Gaia epoch (2016.0 for DR3) to reference observation epoch + init_host_coord_at_obs = init_host_coord.apply_space_motion(Time(self.ref_epoch, format='jd')) + # apply measurement of relative astrometry to get candidate in absolute coordinates at time of observation + init_cand_coord_at_obs = init_host_coord_at_obs.spherical_offsets_by(self.ras[self.ref_epoch_idx] * u.mas, self.decs[self.ref_epoch_idx] * u.mas) + init_cand_coord = SkyCoord(ra=init_cand_coord_at_obs.ra, dec=init_cand_coord_at_obs.dec, distance=Distance(parallax=self.paro*u.mas), pm_ra_cosdec=self.pmrao*u.mas/u.yr, pm_dec=self.pmdeco*u.mas/u.yr, obstime=Time(self.ref_epoch, format='jd')) + # propagate candidate position from epoch of observation to Gaia epoch + tdiff = Time(self.ref_epoch, format='jd').decimalyear-self.gaia_epoch + init_cand_coord_at_gaia = init_cand_coord.apply_space_motion(Time(tdiff+self.gaia_epoch,format='jyear',scale='tcb')) + + # not sure which of these is used anymore + self.ra0 = init_cand_coord_at_gaia.ra.value + self.dec0 = init_cand_coord_at_gaia.dec.value self.pmra0 = self.pmrao # mas/yr self.pmdec0 = self.pmdeco # mas/yr self.par0 = self.paro/10 # mas self.radvel0 = 0 # km/s + # create parameter set for stationary case + # params = ra, dec, pmra, pmdec, par, ra_host, dec_host, pmra_host, pmdec_host, par_host, rv_host + self.stationary_params = [self.ra0, self.dec0, 0, 0, 0, self.rao, self.deco, self.pmrao, self.pmdeco, self.paro, self.radvelo] + # Compute useful chi2 value + self.stationary_loglike = self.loglike(self.stationary_params) + self.stationary_chi2_red = np.log(-2.*self.stationary_loglike)/(len(self.epochs)-self.ndim) + jd_start, jd_end, number = ephem_open() print('[BACKTRACK INFO]: Opened ephemeris file') # if the novas_de405 package is installed this will load the ephemerids file, @@ -529,28 +553,38 @@ def fit(self, dlogz=0.5, npool=4, dynamic=False, nlive=200, mpi_pool=False, resu self.run_median = np.median(samples, axis=0) self.run_quant = np.quantile(samples, [0.1587, 0.8413], axis=0) + # Compute useful chi2 value + self.median_loglike = self.loglike(self.run_median) + self.median_chi2_red = np.log(-2.*self.median_loglike)/(len(self.epochs)-self.ndim) + return self.results def save_results(self, fileprefix='./'): save_dict = {'med': self.run_median, 'quant': self.run_quant, 'results': self.results} target_label = self.target_name.replace(' ','_') file_name = f'{fileprefix}{target_label}_dynestyrun_results.pkl' + print('[BACKTRACK INFO]: Saving results to {}'.format(file_name)) pickle.dump(save_dict, open(file_name, "wb")) def load_results(self, fileprefix: str = './'): target_label = self.target_name.replace(' ', '_') file_name = f'{fileprefix}{target_label}_dynestyrun_results.pkl' + print('[BACKTRACK INFO]: Loading results from {}'.format(file_name)) save_dict = pickle.load(open(file_name, "rb")) self.run_median = save_dict['med'] self.run_quant = save_dict['quant'] self.results = save_dict['results'] + # recompute useful chi2 value + self.median_loglike = self.loglike(self.run_median) + self.median_chi2_red = np.log(-2.*self.median_loglike)/(len(self.epochs)-self.ndim) + def generate_plots( self, days_backward: float = 3.*365., days_forward: float = 3.*365., step_size: float = 10., - ref_epoch: Optional[Tuple[int, int, int]] = None, + # ref_epoch: Optional[Tuple[int, int, int]] = None, plot_radec: bool = False, plot_stationary: bool = False, fileprefix: str = './', @@ -558,12 +592,8 @@ def generate_plots( ) -> Tuple[Figure, Figure, Figure, Figure, Figure]: """ """ - if ref_epoch is None: - mean_epoch = np.floor(np.mean(self.epochs)) - ref_epoch = Time(mean_epoch, format='jd') - - else: - ref_epoch = Time(f'{ref_epoch[0]}-{ref_epoch[1]}-{ref_epoch[2]}T12:00:00', format='isot') + print('[BACKTRACK INFO]: Generating Plots') + ref_epoch = Time(self.ref_epoch, format='jd') fig_track = trackplot( self, @@ -581,5 +611,33 @@ def generate_plots( fig_diag = diagnostic(self, fileprefix=fileprefix, filepost=filepost) fig_prior = plx_prior(self, fileprefix=fileprefix, filepost=filepost) fig_hood = neighborhood(self, fileprefix=fileprefix, filepost=filepost) - + print('[BACKTRACK INFO]: Plots saved to {}'.format(fileprefix)) return fig_track, fig_post, fig_diag, fig_prior, fig_hood + + def generate_stationary_plot( + self, + days_backward: float = 3.*365., + days_forward: float = 3.*365., + step_size: float = 10., + # ref_epoch: Optional[Tuple[int, int, int]] = None, + plot_radec: bool = False, + fileprefix: str = './', + filepost: str = '.pdf', + ) -> Figure: + """ + """ + print('[BACKTRACK INFO]: Generating Stationary plot') + ref_epoch = Time(self.ref_epoch, format='jd') + + fig_track = stationtrackplot( + self, + ref_epoch=ref_epoch.jd, + days_backward=days_backward, + days_forward=days_forward, + step_size=step_size, + plot_radec=plot_radec, + fileprefix=fileprefix, + filepost=filepost + ) + print('[BACKTRACK INFO]: Stationary plot saved to {}'.format(fileprefix)) + return fig_track diff --git a/backtracks/plotting.py b/backtracks/plotting.py index 4a10195..ee35957 100644 --- a/backtracks/plotting.py +++ b/backtracks/plotting.py @@ -135,11 +135,11 @@ def trackplot( # Plot stationary background track at infinity if plot_stationary: - stat_pars = backtracks.run_median.copy() - stat_pars[2:5] = 0.0 + stat_pars = backtracks.stationary_params ra_stat, dec_stat = backtracks.radecdists(plot_epochs_tt, stat_pars) - axs['A'].plot(ra_stat, dec_stat, color="lightgray", label="stationary background", ls='--') + axs['A'].plot(ra_stat, dec_stat, color="lightgray", ls='--', + label="Stationary track, $\chi^2_r={}$".format(round(backtracks.stationary_chi2_red,2))) if plot_radec: axs['B'].plot(plot_times.decimalyear, ra_stat, color="lightgray", ls='--') @@ -193,7 +193,7 @@ def trackplot( ra_bg, dec_bg = backtracks.radecdists(plot_epochs_tt, backtracks.run_median) - axs['A'].plot(ra_bg, dec_bg, color="black", label="best model") + axs['A'].plot(ra_bg, dec_bg, color="black", label="Median track, $\chi^2_r={}$".format(round(backtracks.median_chi2_red,2))) if plot_radec: axs['B'].plot(plot_times.decimalyear, ra_bg, color='black') @@ -223,7 +223,7 @@ def trackplot( axs['A'].errorbar(backtracks.ras, backtracks.decs, yerr=backtracks.decserr, xerr=backtracks.raserr, color="tab:gray", ecolor='tomato', mec='tomato', - label="data", linestyle="none", marker="o", ms=5., mew=1.5) + label="Data", linestyle="none", marker="o", ms=5., mew=1.5) # Plot deltaRA/deltaDec or sep/PA as function of date @@ -341,73 +341,154 @@ def neighborhood(backtracks, fileprefix='./', filepost='.pdf'): return fig -def stationtrackplot(backtracks, ref_epoch, daysback=2600, daysforward=1200, fileprefix='./', filepost='.pdf'): +def stationtrackplot( + backtracks, + ref_epoch, + days_backward=3.*365., + days_forward=3.*365., + step_size=10., + plot_radec=False, + fileprefix='./', + filepost='.pdf' + ): """ """ - fig,axs = plt.subplot_mosaic( + fig, axs = plt.subplot_mosaic( ''' AABB AACC ''', - figsize=(16,8)) + figsize=(16, 8)) - plot_epochs=np.arange(daysforward)+ref_epoch-daysback # 4000 days of epochs to evaluate position at + # Epochs for the background tracks + + plot_epochs = ref_epoch + np.arange(-days_backward, days_forward, step_size) plot_epochs_tt=utc2tt(plot_epochs) + # Create astropy times for observations and background model plot_times = Time(plot_epochs, format='jd') obs_times = Time(backtracks.epochs, format='jd') - post_samples = backtracks.results.samples - best_pars = np.array(backtracks.run_median).T[0] - quant_pars = np.array(backtracks.run_quant).T + # The corr_terms are True where rho is not a nan + + corr_terms =~ np.isnan(backtracks.rho) + + # Plot stationary background track at infinity + + stat_pars = backtracks.stationary_params - corr_terms=~np.isnan(backtracks.rho) # corr_terms is everywhere where rho is not a nan. + ra_stat, dec_stat = backtracks.radecdists(plot_epochs_tt, stat_pars) + axs['A'].plot(ra_stat, dec_stat, color="lightgray", ls='--', + label="Stationary track, $\chi^2_r={}$".format(round(backtracks.stationary_chi2_red,2))) - # plot stationary bg track at infinity (0 parallax) - rasbg,decbg = backtracks.radecdists(plot_epochs_tt, best_pars[0],best_pars[1],0,0,0,best_pars[5],best_pars[6],best_pars[7],best_pars[8],best_pars[9],best_pars[10]) # retrieve coordinates at full range of epochs - # rasbg,decbg = backtracks.radecdists(plot_epochs_tt, backtracks.ra0, backtracks.dec0, 0,0,0) # retrieve coordinates at full range of epochs + if plot_radec: + axs['B'].plot(plot_times.decimalyear, ra_stat, color="lightgray", ls='--') + axs['C'].plot(plot_times.decimalyear, dec_stat, color="lightgray", ls='--') + else: + sep_stat, pa_stat = radec2seppa(ra_stat, dec_stat) + axs['B'].plot(plot_times.decimalyear, sep_stat, color='lightgray', ls='--') + axs['C'].plot(plot_times.decimalyear, pa_stat, color='lightgray', ls='--') + + # Connect data points with best-fit model epochs - axs['B'].plot(plot_times.decimalyear,radec2seppa(rasbg,decbg)[0],color='gray',alpha=1, zorder=3,ls='--') - axs['C'].plot(plot_times.decimalyear,radec2seppa(rasbg,decbg)[1],color='gray',alpha=1, zorder=3,ls='--') + ra_bg_best, dec_bg_best = backtracks.radecdists(backtracks.epochs_tt, stat_pars) - axs['A'].plot(rasbg,decbg,color="gray",zorder=3,label="stationary bg",alpha=1,ls='--') + for i in range(len(backtracks.ras)): + comp_ra = (backtracks.ras[i], ra_bg_best[i]) + comp_dec = (backtracks.decs[i], dec_bg_best[i]) + axs['A'].plot(comp_ra, comp_dec, ls='-', color='tab:gray', lw=1.0) + + # Plot coordinates at observation epochs for best-fit parameters + + axs['A'].plot(ra_bg_best, dec_bg_best, color="tomato", mec='tab:gray', + ms=5., mew=1.5, linestyle="none", marker="o") + + # Plot data points (deltaRA, deltaDEC) + + axs['A'].errorbar(backtracks.ras, backtracks.decs, + yerr=backtracks.decserr, xerr=backtracks.raserr, + color="tab:gray", ecolor='tomato', mec='tomato', + label="Data", linestyle="none", marker="o", ms=5., mew=1.5) + + # Plot deltaRA/deltaDec or sep/PA as function of date + + if plot_radec: + # Plot the deltaRA and deltaDec data points + + axs['B'].errorbar(obs_times.decimalyear, backtracks.ras, + yerr=backtracks.raserr, color="tab:gray", + ecolor='tomato', linestyle="none", + marker='o', ms=5., mew=1.5, mec='tomato') + + axs['C'].errorbar(obs_times.decimalyear, backtracks.decs, + yerr=backtracks.decserr, color="tab:gray", + ecolor='tomato', linestyle="none", + marker='o', ms=5., mew=1.5, mec='tomato') + + else: + # Plot the sep and PA data points + # The error calculation is adopted from orbitize! + + if np.sum(~corr_terms) > 0: + # Plot the sep and PA data points that are not corr_terms + + obs_sep, obs_pa = radec2seppa(backtracks.ras[~corr_terms], backtracks.decs[~corr_terms]) + obs_sep_err = 0.5*backtracks.raserr[~corr_terms] + 0.5*backtracks.decserr[~corr_terms] + obs_pa_err = np.degrees(obs_sep_err/obs_sep) + + axs['B'].errorbar(obs_times[~corr_terms].decimalyear, obs_sep, + yerr=obs_sep_err, color="tab:gray", + ecolor='tomato', linestyle="none", + marker='o', ms=5., mew=1.5, mec='tomato') + + axs['C'].errorbar(obs_times[~corr_terms].decimalyear, obs_pa, + yerr=obs_pa_err, color="tab:gray", + ecolor='tomato', linestyle="none", + marker='o', ms=5., mew=1.5, mec='tomato') + + if np.sum(corr_terms) > 0: + # Plot the sep and PA data points that are corr_terms + + obs_sep, obs_pa = radec2seppa(backtracks.ras[corr_terms], backtracks.decs[corr_terms]) + + obs_sep_err = np.zeros(obs_sep.size) + obs_pa_err = np.zeros(obs_sep.size) + + for i in np.arange(np.sum(corr_terms)): + # Transform the uncertainties from RA/Dec to sep/PA + obs_sep_err[i], obs_pa_err[i], _ = transform_errors( + backtracks.ras[corr_terms][i], backtracks.decs[corr_terms][i], + backtracks.raserr[corr_terms][i], backtracks.decserr[corr_terms][i], + backtracks.rho[corr_terms][i], radec2seppa) + + axs['B'].errorbar(obs_times[corr_terms].decimalyear, obs_sep, + yerr=obs_sep_err, color="tab:gray", + ecolor='tomato', linestyle="none", + marker='o', ms=5., mew=1.5, mec='tomato') + + axs['C'].errorbar(obs_times[corr_terms].decimalyear, obs_pa, + yerr=obs_pa_err, color="tab:gray", + ecolor='tomato', linestyle="none", + marker='o', ms=5., mew=1.5, mec='tomato') - # plot data in deltaRA, deltaDEC plot - axs['A'].errorbar(backtracks.ras[~corr_terms],backtracks.decs[~corr_terms],yerr=backtracks.decserr[~corr_terms],xerr=backtracks.raserr[~corr_terms], - color="tomato",label="data",linestyle="", zorder=5) - axs['A'].errorbar(backtracks.ras[corr_terms],backtracks.decs[corr_terms],yerr=backtracks.decserr[corr_terms],xerr=backtracks.raserr[corr_terms], - color="tomato",linestyle="",marker=".", zorder=5) - # labels - # plt.suptitle("fitted model of background star for YSES-2") axs['A'].invert_xaxis() axs['A'].axis('equal') axs['A'].set_xlabel("Delta RA (mas)") axs['A'].set_ylabel("Delta DEC (mas)") + axs['A'].legend() axs['B'].set_xlabel('Date (year)') axs['C'].set_xlabel('Date (year)') - axs['B'].set_ylabel('Separation (mas)') - axs['C'].set_ylabel('PA (degrees)') - - # plot the datapoints that are not corr_terms in the sep and PA plots, error calculation taken from Orbitize! - sep,pa=radec2seppa(backtracks.ras[~corr_terms],backtracks.decs[~corr_terms]) - sep_err=0.5*backtracks.raserr[~corr_terms]+0.5*backtracks.decserr[~corr_terms] - pa_err=np.degrees(sep_err/sep) - axs['B'].errorbar(obs_times.decimalyear[~corr_terms],sep,yerr=sep_err,color="tomato",linestyle="", zorder=5) - axs['C'].errorbar(obs_times.decimalyear[~corr_terms],pa,yerr=pa_err,color="tomato",linestyle="", zorder=5) - - # plot the corr_terms datapoints with a conversion in the sep and PA plots - sep, pa = radec2seppa(backtracks.ras[corr_terms],backtracks.decs[corr_terms]) - for i in np.arange(np.sum(corr_terms)): - sep_err,pa_err,rho2 = transform_errors(backtracks.ras[corr_terms][i], backtracks.decs[corr_terms][i], - backtracks.raserr[corr_terms][i], backtracks.decserr[corr_terms][i], backtracks.rho[corr_terms][i], radec2seppa) - axs['B'].errorbar(obs_times.decimalyear[corr_terms][i], sep[i], yerr=sep_err, color="tomato", linestyle="", marker='.', zorder=5) - axs['C'].errorbar(obs_times.decimalyear[corr_terms][i], pa[i], yerr=pa_err, color="tomato", linestyle="", marker='.', zorder=5) - axs['A'].legend() + if plot_radec: + axs['B'].set_ylabel('Delta RA (mas)') + axs['C'].set_ylabel('Delta DEC (mas)') + else: + axs['B'].set_ylabel('Separation (mas)') + axs['C'].set_ylabel('PA (degrees)') plt.tight_layout() target_name = backtracks.target_name.replace(' ', '_') - plt.savefig(f"{fileprefix}{target_name}_model_stationary_backtracks"+filepost, dpi=300, bbox_inches='tight') + plt.savefig(f"{fileprefix}{target_name}_stationary_backtracks"+filepost, dpi=300, bbox_inches='tight') return fig