From 4c632def04f4fa286c5a3d93c0982ffe46f1c2ac Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 8 Oct 2024 14:35:59 -0700 Subject: [PATCH 01/13] decimalyear -> julian year --- orbitize/gaia.py | 25 +++++++++++----------- orbitize/hipparcos.py | 50 ++++++++++++++++++++++--------------------- orbitize/sampler.py | 9 +++----- tests/test_gaia.py | 16 +++++++------- 4 files changed, 50 insertions(+), 50 deletions(-) diff --git a/orbitize/gaia.py b/orbitize/gaia.py index af957b7c..d232610d 100644 --- a/orbitize/gaia.py +++ b/orbitize/gaia.py @@ -7,7 +7,7 @@ from astroquery.gaia import Gaia from astropy import units as u import astropy.io.fits as fits -import astropy.time as time +from astropy.time import Time from astropy.io.ascii import read from astropy.coordinates import get_body_barycentric_posvel import numpy.linalg @@ -61,12 +61,12 @@ def __init__(self, gaia_num, hiplogprob, dr="dr2", query=True, gaia_data=None): self.dr = dr if self.dr == "edr3": - self.gaia_epoch = 2016.0 + self.gaia_epoch = Time("J2016", format='jyear_str') elif self.dr == "dr2": - self.gaia_epoch = 2015.5 + self.gaia_epoch = Time("J2015.5", format='jyear_str') else: raise ValueError("`dr` must be either `dr2` or `edr3`") - self.hipparcos_epoch = 1991.25 + self.hipparcos_epoch = Time("J1991.25", format='jyear_str') if query: query = """SELECT @@ -134,12 +134,12 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): """ alpha_H0 = samples[param_idx["alpha0"]] # [deg] - pm_ra = samples[param_idx["pm_ra"]] # [mas/yr] - delta_alpha_from_pm = pm_ra * (self.gaia_epoch - self.hipparcos_epoch) # [mas] + pm_ra = samples[param_idx["pm_ra"]] # [mas/jyr] + delta_alpha_from_pm = pm_ra * (self.gaia_epoch.jyear - self.hipparcos_epoch.jyear) # [mas] delta_H0 = samples[param_idx["delta0"]] # [deg] - pm_dec = samples[param_idx["pm_dec"]] # [mas/yr] - delta_delta_from_pm = pm_dec * (self.gaia_epoch - self.hipparcos_epoch) # [mas] + pm_dec = samples[param_idx["pm_dec"]] # [mas/jyr] + delta_delta_from_pm = pm_dec * (self.gaia_epoch.jyear - self.hipparcos_epoch.jyear) # [mas] # difference in position due to orbital motion between Hipparcos & Gaia epochs alpha_diff_orbit = raoff_model[1, :] - raoff_model[0, :] # [mas] @@ -162,6 +162,7 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): # technically this is an angle so we should wrap it, but the precision # of Hipparcos and Gaia is so good that we'll never have to. alpha_resid = alpha_model - alpha_data + alpha_chi2 = (alpha_resid / alpha_unc) ** 2 delta_model = self.hiplogprob.delta0 + self.mas2deg * ( # [deg] @@ -296,9 +297,9 @@ def __init__(self, hip_id, hiplogprob, gost_filepath, hgca_filepath=None): self.gaia_epoch_dec = entry["epoch_dec_gaia"][0] # read in the GOST file to get the estimated Gaia epochs and scan angles gost_dat = read(gost_filepath, converters={"*": [int, float, bytes]}) - self.gaia_epoch = time.Time( + self.gaia_epoch = Time( gost_dat["ObservationTimeAtGaia[UTC]"] - ).decimalyear # in decimal year + ) # in julian year gaia_scan_theta = np.array(gost_dat["scanAngle[rad]"]) gaia_scan_phi = gaia_scan_theta + np.pi / 2 self.gaia_cos_phi = np.cos(gaia_scan_phi) @@ -443,9 +444,9 @@ def _linear_pm_fit( """ # Sovle y = A * x # construct A matrix - A_pmra = cos_phi * (epochs - epoch_ref_ra) / errs + A_pmra = cos_phi * (epochs.jyear - epoch_ref_ra) / errs A_raoff = cos_phi / errs - A_pmdec = sin_phi * (epochs - epoch_ref_dec) / errs + A_pmdec = sin_phi * (epochs.jyear - epoch_ref_dec) / errs A_decoff = sin_phi / errs A_matrix = np.vstack((A_raoff, A_decoff, A_pmra, A_pmdec)).T diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 72bc4084..0733e6a6 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -1,5 +1,4 @@ import numpy as np -from astropy.io import ascii import pandas as pd import emcee from scipy.stats import norm @@ -17,29 +16,28 @@ class PMPlx_Motion(object): parallax and proper motion model (NO orbital motion is added in this class). Args: - times_mjd (np.array of float): times (in mjd) at which we have absolute astrometric + epochs_mjd (np.array of float): times (in mjd) at which we have absolute astrometric measurements alpha0 (float): measured RA position (in degrees) of the object at alphadec0_epoch (see below). delta0 (float): measured Dec position (in degrees) of the object at alphadec0_epoch (see below). alphadec0_epoch (float): a (fixed) reference time. For stars with Hipparcos data, this - should generally be 1991.25, but you can define it however you want. Absolute + should generally be J1991.25, but you can define it however you want. Absolute astrometric data (passed in via an orbitize! data table) should be defined as offsets from the reported position of the object at this epoch (with propagated uncertainties). For example, if you have two absolute astrometric measurements - taken with GRAVITY, as well as a Hipparcos-derived position (at epoch 1991.25), - alphadec0_epoch should be 1991.25, and you should pass in absolute astrometry + taken with GRAVITY, as well as a Hipparcos-derived position (at epoch J1991.25), + alphadec0_epoch should be J1991.25, and you should pass in absolute astrometry in terms of mas *offset* from the Hipparcos catalog position, with propagated errors of your measurement and the Hipparcos measurement. """ - def __init__(self, epochs_mjd, alpha0, delta0, alphadec0_epoch=1991.25): + def __init__(self, epochs_mjd, alpha0, delta0, alphadec0_epoch=Time("J1991.25", format="jyear_str").jyear): self.epochs_mjd = epochs_mjd self.alphadec0_epoch = alphadec0_epoch self.alpha0 = alpha0 self.delta0 = delta0 epochs = Time(epochs_mjd, format="mjd") - self.epochs = epochs.decimalyear # compute Earth XYZ position in barycentric coordinates bary_pos, _ = get_body_barycentric_posvel("earth", epochs) @@ -61,7 +59,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): indices in an array of fitting parameters (generally set to System.basis.param_idx). epochs: if None, use self.epochs for astrometric predictions. Otherwise, - use this array passed in [in decimalyear]. + use this array passed in [in mjd]. Returns: tuple of: @@ -78,14 +76,14 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): delta_H0 = samples[param_idx["delta0"]] if epochs is None: - epochs = self.epochs + epochs = Time(self.epochs_mjd, format='mjd') X = self.X Y = self.Y Z = self.Z else: # compute Earth XYZ position in barycentric coordinates bary_pos, _ = get_body_barycentric_posvel( - "earth", Time(epochs, format="decimalyear") + "earth", Time(epochs, format="mjd") ) X = bary_pos.x.value # [au] Y = bary_pos.y.value # [au] @@ -105,7 +103,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): X[i] * np.sin(np.radians(self.alpha0)) - Y[i] * np.cos(np.radians(self.alpha0)) ) - + (epochs[i] - self.alphadec0_epoch) * pm_ra + + (epochs[i].jyear - self.alphadec0_epoch) * pm_ra ) delta_C_array[i] = ( delta_H0 @@ -119,7 +117,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): * np.sin(np.radians(self.delta0)) - Z[i] * np.cos(np.radians(self.delta0)) ) - + (epochs[i] - self.alphadec0_epoch) * pm_dec + + (epochs[i].jyear - self.alphadec0_epoch) * pm_dec ) return alpha_C_st_array, delta_C_array @@ -137,9 +135,9 @@ class HipparcosLogProb(object): are described here for completeness. See Nielsen+ 2020 for more detail. - alpha0: RA offset from the reported Hipparcos position at a particular - epoch (usually 1991.25) [mas] + epoch (usually J1991.25) [mas] - delta0: Dec offset from the reported Hipparcos position at a particular - epoch (usually 1991.25) [mas] + epoch (usually J1991.25) [mas] - pm_ra: RA proper motion [mas/yr] - pm_dec: Dec proper motion [mas/yr] - plx: parallax [mas] @@ -157,7 +155,7 @@ class HipparcosLogProb(object): zeros in the prefix if number is <100,000. (i.e. 27321 should be passed in as '027321'). num_secondary_bodies (int): number of companions in the system - alphadec0_epoch (float): epoch (in decimal year) that the fitting + alphadec0_epoch (float): epoch (in Julian decimal year) that the fitting parameters alpha0 and delta0 are defined relative to (see above). renormalize_errors (bool): if True, normalize the scan errors to get chisq_red = 1, following Nielsen+ 2020 (eq 10). In general, this @@ -173,7 +171,7 @@ def __init__( path_to_iad_file, hip_num, num_secondary_bodies, - alphadec0_epoch=1991.25, + alphadec0_epoch=Time("J1991.25", format="jyear_str").jyear, renormalize_errors=False, ): self.path_to_iad_file = path_to_iad_file @@ -284,7 +282,11 @@ def __init__( n_lines = len(iad) - times = iad[1] + 1991.25 + + self.epochs = Time( + iad[1] * 365.25 + Time("J1991.25",format="jyear_str").jd, format='jd' + ) + self.cos_phi = iad[3] # scan direction self.sin_phi = iad[4] self.R = iad[5] # abscissa residual [mas] @@ -295,18 +297,18 @@ def __init__( if n_lines - len(good_scans) > 0: print("{} Hipparcos scans rejected.".format(n_lines - len(good_scans))) - times = times[good_scans] + self.epochs = self.epochs[good_scans] self.cos_phi = self.cos_phi[good_scans] self.sin_phi = self.sin_phi[good_scans] self.R = self.R[good_scans] self.eps = self.eps[good_scans] + self.epochs_mjd = self.epochs.mjd + # if the star has a type 1 (stochastic) solution, we need to undo the addition of a jitter term in quadrature self.eps = np.sqrt(self.eps**2 - self.var**2) - epochs = Time(times, format="decimalyear") - self.epochs = epochs.decimalyear - self.epochs_mjd = epochs.mjd + self.hipparcos_plxpm_predictor = PMPlx_Motion( self.epochs_mjd, @@ -316,7 +318,7 @@ def __init__( ) if self.renormalize_errors: - D = len(epochs) - 6 + D = len(self.epochs) - 6 G = f2 f = (G * np.sqrt(2 / (9 * D)) + 1 - (2 / (9 * D))) ** (3 / 2) @@ -333,7 +335,7 @@ def __init__( self.hipparcos_plxpm_predictor.X * np.sin(np.radians(self.alpha0)) - self.hipparcos_plxpm_predictor.Y * np.cos(np.radians(self.alpha0)) ) - + (self.epochs - 1991.25) * self.pm_ra0 + + (self.epochs.jyear - self.alphadec0_epoch) * self.pm_ra0 ) changein_delta = ( @@ -347,7 +349,7 @@ def __init__( * np.sin(np.radians(self.delta0)) - self.hipparcos_plxpm_predictor.Z * np.cos(np.radians(self.delta0)) ) - + (self.epochs - 1991.25) * self.pm_dec0 + + (self.epochs.jyear - self.alphadec0_epoch) * self.pm_dec0 ) # compute abcissa point (Nielsen+ Eq 3) diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 5ee50f21..e9ad25e2 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -126,12 +126,9 @@ def _logl(self, params): ) if self.system.gaia is not None: - gaiahip_epochs = Time( - np.append( - self.system.gaia.hipparcos_epoch, self.system.gaia.gaia_epoch - ), - format="decimalyear", - ).mjd + gaiahip_epochs = np.append( + self.system.gaia.hipparcos_epoch.mjd, self.system.gaia.gaia_epoch.mjd + ) # compute Ra/Dec predictions at the Gaia epoch raoff_model, deoff_model, _ = self.system.compute_all_orbits( diff --git a/tests/test_gaia.py b/tests/test_gaia.py index 70dbce85..ebc0e840 100644 --- a/tests/test_gaia.py +++ b/tests/test_gaia.py @@ -156,11 +156,11 @@ def test_orbit_calculation(): myGaia.ra = myHip.alpha0 + ( myGaia.mas2deg * pm_a - * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) + * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear) / np.cos(np.radians(myHip.delta0)) ) myGaia.dec = myHip.delta0 + ( - myGaia.mas2deg * pm_d * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) + myGaia.mas2deg * pm_d * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear) ) test_samples = [sma, ecc, inc, aop, pan, tau, plx, m0, m1, a0, d0, pm_a, pm_d] @@ -194,7 +194,7 @@ def test_orbit_calculation(): myGaia.dec = myHip.delta0 + 1 sma = 2 * (myGaia.dec - myHip.delta0) * deg2arcsec * (plx * mas2arcsec) # [au] - per = 2 * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) # [yr] + per = 2 * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear) # [yr] mtot = sma**3 / per**2 test_samples[param_idx["sma1"]] = sma @@ -203,7 +203,7 @@ def test_orbit_calculation(): # passes through peri (+sma decl for e=0 orbits) at Hipparcos epoch # -> @ Gaia epoch, primary should be at +sma decl - tau = basis.tp_to_tau(myGaia.hipparcos_epoch, 58849, per) + tau = basis.tp_to_tau(myGaia.hipparcos_epoch.jyear, 58849, per) test_samples[param_idx["tau1"]] = tau # choose sma and mass so that Hipparcos/Gaia difference is only due to orbit. @@ -308,8 +308,8 @@ def test_nointernet(): if __name__ == "__main__": test_nointernet() - # test_dr2_edr3() - # test_system_setup() - # test_valueerror() - # test_orbit_calculation() + test_dr2_edr3() + test_system_setup() + test_valueerror() + test_orbit_calculation() test_hgca() From 3feb266c43753964ee3bba8116ad34f9e136eb82 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 15 Jun 2026 09:49:15 -0700 Subject: [PATCH 02/13] fix unit tests --- tests/test_abs_astrometry.py | 9 ++++----- tests/test_hipparcos.py | 6 +++--- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/tests/test_abs_astrometry.py b/tests/test_abs_astrometry.py index ca6c3a8a..88626890 100644 --- a/tests/test_abs_astrometry.py +++ b/tests/test_abs_astrometry.py @@ -88,10 +88,9 @@ def test_arbitrary_abs_astrom(): path_to_iad_file, hip_num, num_secondary_bodies ) - epochs_astropy = Time( - np.array([0, 0.5, 1.0]) + testHiPIAD.alphadec0_epoch, format="decimalyear" - ) - epochs = epochs_astropy.mjd + epochs = Time( + testHiPIAD.alphadec0_epoch, format="jyear" + ).mjd + np.array([0, 365.25/2, 365.25]) ra_model = np.zeros(epochs.shape) dec_model = np.zeros(epochs.shape) @@ -190,5 +189,5 @@ def test_arbitrary_abs_astrom(): if __name__ == "__main__": - test_1planet() + # test_1planet() test_arbitrary_abs_astrom() diff --git a/tests/test_hipparcos.py b/tests/test_hipparcos.py index 08299835..5c401a47 100644 --- a/tests/test_hipparcos.py +++ b/tests/test_hipparcos.py @@ -148,7 +148,7 @@ def test_dvd_vs_2021catalog(): np.isclose(new_iadHipLogProb.sin_phi, old_iadHipLogProb.sin_phi, atol=1e-2) ) assert np.all( - np.isclose(new_iadHipLogProb.epochs, old_iadHipLogProb.epochs, atol=1e-2) + np.isclose(new_iadHipLogProb.epochs.value, old_iadHipLogProb.epochs.value, atol=1e-2) ) @@ -268,5 +268,5 @@ def test_save_load_2021(): # test_save_load_dvd() # test_save_load_2021() # test_hipparcos_api() - test_iad_refitting() - # test_dvd_vs_2021catalog() + # test_iad_refitting() + test_dvd_vs_2021catalog() From d341b3cc0d30a3b7e3915d8b6f490f554e819db3 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 15 Jun 2026 10:01:18 -0700 Subject: [PATCH 03/13] fix proper motion plot --- orbitize/plot.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/orbitize/plot.py b/orbitize/plot.py index 7f248695..f8c6f6f7 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -1673,7 +1673,7 @@ def plot_propermotion( axs[0].set_ylabel(r"$\mu_\alpha^*$ [mas/yr]") axs[0].errorbar( - np.nanmedian(system.gaia.hipparcos_epoch), + np.nanmedian(system.gaia.hipparcos_epoch.jyear), system.gaia.hip_pm[0], yerr=system.gaia.hip_pm_err[0], zorder=30, @@ -1683,10 +1683,10 @@ def plot_propermotion( ) hgca_epoch = ( - system.gaia.gaia_epoch_ra + np.nanmedian(system.gaia.hipparcos_epoch) + system.gaia.gaia_epoch_ra + np.nanmedian(system.gaia.hipparcos_epoch.jyear) ) / 2 hgca_epoch_err = ( - system.gaia.gaia_epoch_ra - np.nanmedian(system.gaia.hipparcos_epoch) + system.gaia.gaia_epoch_ra - np.nanmedian(system.gaia.hipparcos_epoch.jyear) ) / 2 axs[0].errorbar( @@ -1714,7 +1714,7 @@ def plot_propermotion( axs[1].yaxis.set_major_formatter(FormatStrFormatter("%.1f")) axs[1].errorbar( - np.nanmedian(system.gaia.hipparcos_epoch), + np.nanmedian(system.gaia.hipparcos_epoch.jyear), system.gaia.hip_pm[1], yerr=system.gaia.hip_pm_err[1], zorder=30, From 8842f4713750d9846bfc6fa87e45a4bd37a25f7d Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 28 Jul 2026 13:15:16 -0400 Subject: [PATCH 04/13] use float for alphadec0_epoch input --- orbitize/hipparcos.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index b1fbab54..94908548 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -20,7 +20,7 @@ class PMPlx_Motion(object): measurements alpha0 (float): measured RA position (in degrees) of the object at alphadec0_epoch (see below). delta0 (float): measured Dec position (in degrees) of the object at alphadec0_epoch (see below). - alphadec0_epoch (float): a (fixed) reference time. For stars with Hipparcos data, this + alphadec0_epoch (float): a (fixed) reference time (in jyear). For stars with Hipparcos data, this should generally be J1991.25, but you can define it however you want. Absolute astrometric data (passed in via an orbitize! data table) should be defined as offsets from the reported position of the object at this epoch (with propagated @@ -31,7 +31,7 @@ class PMPlx_Motion(object): errors of your measurement and the Hipparcos measurement. """ - def __init__(self, epochs_mjd, alpha0, delta0, alphadec0_epoch=Time("J1991.25", format="jyear_str").jyear): + def __init__(self, epochs_mjd, alpha0, delta0, alphadec0_epoch=1991.25): self.epochs_mjd = epochs_mjd self.alphadec0_epoch = alphadec0_epoch self.alpha0 = alpha0 From d35404ead663c36a537361aafdd225d80847318e Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 28 Jul 2026 13:53:13 -0400 Subject: [PATCH 05/13] draft e2e hgca test --- orbitize/example_data/H024332.d | 94 +++++++++++++++++++ orbitize/example_data/HD33632Ab.csv | 4 + ....3_806543_2026-07-28-17-46-26_HD_33632.csv | 76 +++++++++++++++ tests/end-to-end-tests/hd_33632Ab.py | 57 +++++++++++ 4 files changed, 231 insertions(+) create mode 100644 orbitize/example_data/H024332.d create mode 100644 orbitize/example_data/HD33632Ab.csv create mode 100644 orbitize/example_data/gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv create mode 100644 tests/end-to-end-tests/hd_33632Ab.py diff --git a/orbitize/example_data/H024332.d b/orbitize/example_data/H024332.d new file mode 100644 index 00000000..d7d1682f --- /dev/null +++ b/orbitize/example_data/H024332.d @@ -0,0 +1,94 @@ +# This file contains residual records, extracted from the Hipparcos 2 +# Interactive Data Access Tool (2014). For more information, see: +# https://www.cosmos.esa.int/web/hipparcos/interactive-data-access +# https://www.cosmos.esa.int/web/hipparcos/catalogues +# +# HIP MCE NRES NC isol_n SCE F2 F1 +# 24332 24273 81 1 5 0 1.17 0 +# Hp B-V VarAnn NOB NR +# 6.5960 0.542 0 81 0 +# RAdeg DEdeg Plx pm_RA pm_DE e_RA e_DE e_Plx e_pmRA e_pmDE dpmRA dpmDE e_dpmRA e_dpmDE ddpmRA ddpmDE e_ddpmRA e_ddpmDE upsRA upsDE e_upsRA e_upsDE var +# 78.32313339 37.33763724 38.29 -145.00 -135.14 0.51 0.26 0.55 0.50 0.27 --- --- --- --- --- --- --- --- --- --- --- --- --- +# +# IORB EPOCH PARF CPSI SPSI RES SRES + 191 -1.1740 0.2376 -0.1869 0.9824 0.63 1.81 + 191 -1.1740 0.2375 -0.1869 0.9824 0.70 2.22 + 192 -1.1727 0.2280 -0.1750 0.9846 1.29 1.62 + 192 -1.1727 0.2289 -0.1753 0.9845 0.24 1.94 + 192 -1.1727 0.2252 -0.1706 0.9853 -1.62 1.60 + 192 -1.1727 0.2241 -0.1690 0.9856 3.60 2.00 + 192 -1.1727 0.2220 -0.1662 0.9861 2.60 1.48 + 216 -1.1435 -0.4636 0.5895 0.8078 -1.74 1.73 + 216 -1.1435 -0.4638 0.5896 0.8077 0.24 1.27 + 216 -1.1435 -0.4680 0.5943 0.8043 1.42 1.65 + 290 -1.0536 0.6767 -0.7471 0.6647 -0.89 1.73 + 290 -1.0536 0.6758 -0.7463 0.6656 -4.63 1.81 + 290 -1.0536 0.6776 -0.7479 0.6638 -3.15 1.70 + 290 -1.0536 0.6771 -0.7474 0.6644 0.82 1.96 + 290 -1.0536 0.6782 -0.7485 0.6632 -0.64 1.25 + 371 -0.9551 -0.3388 0.1595 0.9872 -1.68 1.61 + 371 -0.9551 -0.3387 0.1594 0.9872 -1.85 1.40 + 371 -0.9551 -0.3344 0.1540 0.9881 -0.18 1.43 + 372 -0.9539 -0.3252 0.1405 0.9901 -0.38 1.71 + 372 -0.9539 -0.3242 0.1393 0.9902 0.60 1.49 + 372 -0.9539 -0.3207 0.1343 0.9909 -1.39 1.82 + 372 -0.9539 -0.3187 0.1318 0.9913 -0.15 1.74 + 372 -0.9539 -0.3158 0.1277 0.9918 -1.45 1.96 + 372 -0.9539 -0.3145 0.1262 0.9920 -1.43 1.90 + 372 -0.9539 -0.3092 0.1194 0.9928 -0.59 1.91 + 373 -0.9527 -0.3025 0.1093 0.9940 -0.35 1.97 + 373 -0.9527 -0.3013 0.1078 0.9942 -1.47 1.92 + 373 -0.9527 -0.2972 0.1018 0.9948 0.46 1.80 + 385 -0.9381 -0.0032 -0.3088 0.9511 2.46 2.20 + 385 -0.9381 -0.0017 -0.3108 0.9505 0.89 2.34 + 386 -0.9371 0.0058 -0.3229 0.9464 -1.15 2.06 + 386 -0.9371 0.0110 -0.3300 0.9440 0.99 1.85 + 387 -0.9356 0.0269 -0.3539 0.9353 0.37 1.58 + 387 -0.9356 0.0281 -0.3558 0.9346 0.31 1.46 + 387 -0.9356 0.0368 -0.3662 0.9305 1.17 1.65 + 387 -0.9356 0.0371 -0.3676 0.9300 -1.05 1.42 + 677 -0.5832 -0.6617 -0.6426 -0.7662 -0.41 1.91 + 677 -0.5832 -0.6609 -0.6419 -0.7668 -3.56 2.10 + 745 -0.5005 0.6159 0.7602 -0.6497 -1.31 1.25 + 745 -0.5005 0.6152 0.7594 -0.6506 -1.75 1.26 + 745 -0.5005 0.6140 0.7584 -0.6518 -3.36 1.59 + 745 -0.5005 0.6112 0.7559 -0.6547 3.84 1.80 + 783 -0.4543 -0.4032 -0.2511 -0.9680 -0.80 2.98 + 783 -0.4543 -0.4077 -0.2567 -0.9665 -3.22 2.88 + 783 -0.4543 -0.4071 -0.2562 -0.9666 2.74 2.52 + 1071 -0.1044 0.6198 -0.6600 0.7513 0.71 1.83 + 1071 -0.1044 0.6177 -0.6578 0.7532 0.61 1.35 + 1071 -0.1044 0.6183 -0.6584 0.7526 0.60 1.12 + 1128 -0.0351 -0.6614 0.5774 0.8165 -0.30 1.31 + 1128 -0.0351 -0.6621 0.5780 0.8160 0.85 1.28 + 1128 -0.0351 -0.6598 0.5755 0.8178 -3.85 1.85 + 1128 -0.0351 -0.6600 0.5758 0.8176 2.53 1.76 + 1176 0.0232 0.4723 -0.7162 0.6979 0.17 3.08 + 1176 0.0232 0.4724 -0.7163 0.6977 -4.10 3.56 + 1176 0.0232 0.4754 -0.7194 0.6946 1.57 3.41 + 1176 0.0232 0.4761 -0.7195 0.6944 -4.89 3.12 + 1457 0.3647 -0.4952 -0.6027 -0.7980 1.73 3.32 + 1457 0.3647 -0.4940 -0.6016 -0.7988 -1.52 3.16 + 1457 0.3647 -0.4921 -0.5994 -0.8004 -2.06 3.26 + 1457 0.3647 -0.4910 -0.5985 -0.8011 -0.85 3.08 + 1458 0.3658 -0.4904 -0.5940 -0.8045 0.06 3.15 + 1508 0.4266 0.6755 0.7105 -0.7037 0.57 3.49 + 1854 0.8469 0.4172 -0.4309 0.9024 1.16 2.49 + 1854 0.8469 0.4162 -0.4296 0.9030 8.28 2.66 + 1854 0.8469 0.4114 -0.4238 0.9058 1.80 2.28 + 1855 0.8481 0.4098 -0.4218 0.9067 -2.22 1.61 + 1855 0.8481 0.4096 -0.4220 0.9066 -1.14 1.64 + 1892 0.8930 -0.6099 0.6390 0.7692 -2.99 1.25 + 1892 0.8930 -0.6121 0.6410 0.7675 0.48 1.96 + 1892 0.8930 -0.6108 0.6397 0.7686 -1.69 2.03 + 1957 0.9719 0.6484 -0.7590 0.6511 0.47 1.56 + 1958 0.9732 0.6487 -0.7614 0.6482 -1.34 1.27 + 2245 1.3217 -0.0755 -0.2542 -0.9671 3.35 2.64 + 2246 1.3228 -0.0634 -0.2349 -0.9720 2.06 2.23 + 2246 1.3228 -0.0631 -0.2348 -0.9720 0.27 2.16 + 2246 1.3228 -0.0588 -0.2292 -0.9734 1.09 1.51 + 2246 1.3228 -0.0586 -0.2288 -0.9735 -2.37 1.40 + 2246 1.3228 -0.0540 -0.2228 -0.9749 -1.92 1.32 + 2246 1.3228 -0.0535 -0.2220 -0.9751 -1.69 1.34 + 2717 1.8945 0.6215 -0.6608 0.7505 -1.56 2.36 + 2717 1.8945 0.6208 -0.6601 0.7512 -1.05 1.36 diff --git a/orbitize/example_data/HD33632Ab.csv b/orbitize/example_data/HD33632Ab.csv new file mode 100644 index 00000000..3e3fae57 --- /dev/null +++ b/orbitize/example_data/HD33632Ab.csv @@ -0,0 +1,4 @@ +epoch,object,raoff,raoff_err,decoff,decoff_err,inst +58409.0,1,-0.761,0.005,-0.176,0.004,CHARIS +58423.0,1,-0.753,0.005,-0.178,0.005,NIRC2 +59092.0,1,-0.740,0.005,-0.095,0.003,CHARIS \ No newline at end of file diff --git a/orbitize/example_data/gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv b/orbitize/example_data/gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv new file mode 100644 index 00000000..37502d39 --- /dev/null +++ b/orbitize/example_data/gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv @@ -0,0 +1,76 @@ +Target, ra[rad], dec[rad], ra[h:m:s], dec[d:m:s], ObservationTimeAtGaia[UTC], CcdRow[1-7], zetaFieldAngle[rad], scanAngle[rad], Fov[FovP=preceding/FovF=following], parallaxFactorAlongScan, parallaxFactorAcrossScan, ObservationTimeAtBarycentre[BarycentricJulianDateInTCB] +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2014-09-07T20:02:55.945,1,0.004261441385324073,-2.419338622854937,FoVP,-0.7139807546396543,0.721507425468232,2456908.3357978356 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2014-09-07T21:49:30.161,5,-9.544894992248997E-4,-2.4202204741216615,FoVF,-0.7136305588529797,0.7219838864565272,2456908.409812096 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2014-10-11T19:00:12.343,4,0.001844990104153869,-2.6989244787331725,FoVF,-0.5481129956895733,0.715562776608969,2456942.295425042 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-02-26T06:30:37.544,4,-0.0010163849352416324,-0.7842990960205769,FoVP,0.6685663841239335,0.7038981947222167,2457079.7735618614 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-02-26T08:17:11.748,2,0.0043439215473310136,-0.7831533431631229,FoVF,0.6677295026077388,0.7051080241324862,2457079.8475617887 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-03-28T07:55:49.189,1,0.005457876467591252,0.588823806199045,FoVF,-0.6552996516448154,0.7077481533077975,2457109.8298096643 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-03-28T12:09:29.205,7,-0.005608569703553211,0.5823156230017611,FoVP,-0.6507356749560859,0.7108871281275355,2457110.005950688 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-04-18T06:12:04.293,7,-0.006533791344582265,-0.7680639806381448,FoVP,0.3898713902485748,0.7177846905785434,2457130.7559392946 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-04-18T07:58:38.435,6,-0.003103768767807712,-0.7721436609584973,FoVF,0.3921103718707833,0.7157601708853536,2457130.829939783 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-08-31T06:43:00.914,1,0.004865564971320537,-2.4058116280271005,FoVP,-0.6878085837994734,0.7220322118950759,2457265.7795853484 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-08-31T08:29:35.136,5,-1.3801332074555073E-4,-2.407470299670282,FoVF,-0.6870095371493553,0.723135527289165,2457265.8535995437 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-10-21T06:18:54.018,1,0.004471299965674476,-2.8968391351998934,FoVP,-0.40082547260616475,0.7152650870733186,2457316.7675270303 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-10-21T08:05:28.195,4,5.314581390129491E-4,-2.893116607845149,FoVF,-0.40346347720117237,0.7129338131237428,2457316.84153943 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T02:22:49.398,7,-0.006579356942279041,0.236145041021683,FoVP,-0.03834907693935221,0.694428649438883,2457411.6042917464 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T08:23:03.958,6,-0.0038725920968205717,0.2546928790991775,FoVP,-0.052999549896909234,0.6964123003569312,2457411.85444359 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T10:09:38.234,6,-0.0031575038856307617,0.2601916339504578,FoVF,-0.057358402889244674,0.6969374272071714,2457411.9284464642 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T14:23:18.360,4,-0.0016147811680190876,0.27329939400961767,FoVP,-0.06777602079582802,0.6980729604950042,2457412.1045935275 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T16:09:52.595,5,-0.0010341575185062362,0.27881361290059464,FoVF,-0.07216954405264611,0.6985016850502285,2457412.178595899 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T20:23:32.633,3,1.8556436156794566E-4,0.2919545529356563,FoVP,-0.08266471202294118,0.699406370240556,2457412.3547418835 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T22:10:06.835,4,6.292818786538373E-4,0.29748114711269064,FoVF,-0.08708867642283567,0.6997375811664404,2457412.4287438514 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T04:10:20.983,4,0.0018248394530124305,0.3161833991809132,FoVF,-0.10210103063515528,0.7006418920781364,2457412.6788906604 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T08:24:00.906,2,0.0023818943010621763,0.32936757429373015,FoVP,-0.1127198448698961,0.7010784128387068,2457412.8550352077 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T10:10:35.068,3,0.0025450006018311066,0.33490894028439583,FoVF,-0.1171911823082803,0.7012122611944958,2457412.9290366755 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T14:24:14.967,2,0.00276305700835884,0.3481023527273083,FoVP,-0.12785522427154672,0.7014127131565167,2457413.1051808777 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T16:10:49.123,3,0.0027827504038084585,0.3536457808391573,FoVF,-0.13234316126613765,0.7014472524242198,2457413.179182254 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T20:24:29.019,2,0.002657189307499812,0.3668396702564774,FoVP,-0.1430409223748212,0.7014112029115603,2457413.355326366 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T22:11:03.178,3,0.0025316119796161395,0.37238140551992394,FoVF,-0.14754057165499984,0.7013463982839263,2457413.429327759 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T02:24:43.093,2,0.002058215448561157,0.3855666646156358,FoVP,-0.15826032082384694,0.701074122304188,2457413.605472041 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T04:11:17.266,4,0.0017856834223902193,0.39110281243277567,FoVF,-0.16276672159993388,0.7009102418277547,2457413.67947356 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T08:24:57.222,3,9.606596044898878E-4,0.40427001689142494,FoVP,-0.17349660196053956,0.7004027467082757,2457413.855618272 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T10:11:31.418,4,5.396729360711204E-4,0.40979655614814264,FoVF,-0.1780047614673837,0.7001403707809352,2457413.9296200266 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T14:25:11.439,4,-6.403204596819836E-4,0.42293599705504864,FoVP,-0.18873289084544134,0.6993994168962655,2457414.1057654307 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T16:11:45.665,5,-0.0012110681643103644,0.42844879422358995,FoVF,-0.19323782756796007,0.6990394454062332,2457414.179767531 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T20:25:25.773,5,-0.002748902097942923,0.4415505126895518,FoVP,-0.20395240322629904,0.6980675596501031,2457414.3559138877 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T22:12:00.041,7,-0.003470514929325013,0.4470453367865934,FoVF,-0.20844919124041658,0.6976112166242967,2457414.4299164424 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-27T02:25:40.259,6,-0.005368569191757951,0.46009916042917925,FoVP,-0.2191385960607972,0.6964117005468431,2457414.6060640095 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-03-07T19:39:51.913,3,0.0012121922309560186,-0.8542754870348369,FoVP,0.7041869581442349,0.7060162146074431,2457455.320674432 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-03-07T21:26:26.122,1,0.006665552464438079,-0.8542565304598202,FoVF,0.7040186443066007,0.7063161809805514,2457455.3946742327 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-12T14:55:55.710,1,0.005942206012869388,0.3628258985486866,FoVF,-0.48780566568372496,0.7085666572351743,2457491.1201509284 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-12T19:09:35.826,4,-7.294599410350257E-4,0.35199461456885456,FoVP,-0.4799952817800107,0.7120836520516343,2457491.296295099 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-12T20:56:10.002,7,-0.0034748396872719325,0.34740972259971487,FoVF,-0.47668869637983835,0.7135501447550359,2457491.3702956517 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T13:17:49.888,7,-0.00639600613720418,-0.5340341406546021,FoVP,0.1506039943362212,0.7184719412369589,2457504.051040776 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T15:04:24.008,7,-0.004297615005516056,-0.5389084680464333,FoVF,0.15350936115572406,0.7170256827047164,2457504.1250416706 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T19:18:03.957,3,7.591508228285377E-4,-0.5504829060947667,FoVP,0.16036430913827335,0.7135379064159008,2457504.3011863227 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T21:04:38.197,3,0.0029092713056413024,-0.5553339189506192,FoVF,0.1632188447704601,0.7120539222954,2457504.375188631 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-26T16:07:23.284,1,0.006574733665764367,-3.118795104425183,FoVF,0.10682475654855952,0.7218811268715385,2457596.168649828 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-26T22:07:37.867,2,0.005180252835275456,-3.1363485488714855,FoVF,0.11848620507860806,0.7229014193780243,2457596.4188352157 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T02:21:18.090,1,0.004362678166136526,3.134463758074025,FoVP,0.12674718221109432,0.7234981166164549,2457596.5950069344 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T04:07:52.375,2,0.004060258859464345,3.1292629484774115,FoVF,0.13022930699956975,0.7237185581011203,2457596.6690198025 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T08:21:32.554,1,0.0034392270401076255,3.116877979780593,FoVP,0.13854387784111935,0.724170859952826,2457596.8451910713 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T10:08:06.823,3,0.003220174398940672,3.111672752562402,FoVF,0.14204752721561414,0.7243303340855934,2457596.919203779 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T14:21:46.971,2,0.0027994552690540557,3.0992787745764105,FoVP,0.1504109179466158,0.7246368681123053,2457597.095374736 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T16:08:21.229,3,0.002665340738750663,3.0940704121321847,FoVF,0.1539340462749031,0.7247348503782834,2457597.169387342 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T20:22:01.358,2,0.0024486425767270693,3.081670569837051,FoVP,0.16234129836396347,0.7248944717745328,2457597.3455581274 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T22:08:35.611,3,0.0024010093550536505,3.0764604318182083,FoVF,0.16588178662740793,0.7249305255718157,2457597.4195706914 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T02:22:15.733,2,0.002391968504679075,3.0640580527793175,FoVP,0.17432777634488222,0.7249423130042384,2457597.5957414494 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T04:08:49.985,3,0.0024323280144230033,3.058847575269989,FoVF,0.17788344036041448,0.724916102015573,2457597.669754032 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T08:22:30.113,2,0.0026345009197534094,3.046446167753184,FoVP,0.18636289987487514,0.7247793795003656,2457597.8459249088 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T10:09:04.371,3,0.002764329210817298,3.0412368621777817,FoVF,0.1899314997313748,0.7246906732262365,2457597.9199375724 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T14:22:44.517,2,0.003181183779839112,3.0288401122270505,FoVP,0.19843904058849557,0.7244050230103958,2457598.0961087164 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T16:09:18.786,3,0.0034019178732009876,3.0236335639869374,FoVF,0.20201829092939788,0.7242537020485526,2457598.170121524 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T20:22:58.964,1,0.004036824717549616,3.0112453318580044,FoVP,0.2105484294938558,0.7238189768308214,2457598.346293086 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T22:09:33.248,2,0.004349858903315388,3.006043198698712,FoVF,0.21413601060055207,0.7236050378553299,2457598.4203061014 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-29T04:09:47.778,1,0.005612764589441238,2.9884715247663682,FoVF,0.22627676495516516,0.7227449326322577,2457598.6704915203 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-09-10T07:53:21.077,3,2.5697815122738374E-4,-2.4290032247557622,FoVP,-0.7192512483154985,0.7205751339300248,2457641.8294886546 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-17T19:13:04.480,2,0.0019670414371021298,2.3267191540506578,FoVP,0.4528101981918288,0.7141797978884794,2457679.304956016 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-17T20:59:38.636,2,0.0049775136690584825,2.3312166991054286,FoVF,0.44890664506950456,0.7158274316724391,2457679.378968362 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-28T09:06:29.801,3,0.003469437222110181,3.065978956067008,FoVF,-0.180726009925592,0.712925193808823,2457689.884508299 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-28T13:20:09.807,5,-0.002756221733663789,3.0773433937782557,FoVP,-0.18931149002963266,0.7086098902236475,2457690.060677803 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-01-28T21:19:07.703,5,-0.0022362656409328804,-0.26838445629601826,FoVP,0.29847350987287763,0.6975925921145726,2457782.3930357434 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-01-28T23:05:41.920,4,0.0011818264198154163,-0.2640267638820214,FoVF,0.2953490008385252,0.6997892657822943,2457782.4670374626 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-02-12T04:47:43.319,4,-2.97168219826661E-4,0.7062620273648472,FoVP,-0.5481689981488462,0.7034066986201087,2457796.7034035507 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-02-12T06:34:17.557,7,-0.0044351025601507554,0.709904993482966,FoVF,-0.5513672776360866,0.7016248404356326,2457796.7774044797 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-03-18T10:12:07.737,6,-0.004262440417520913,-0.891777649505845,FoVP,0.697256372367284,0.7091677863936382,2457830.925425397 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-03-18T11:58:41.940,4,0.0010442329234670706,-0.8928944307050657,FoVF,0.6977166434167678,0.7085186421825921,2457830.9994251765 diff --git a/tests/end-to-end-tests/hd_33632Ab.py b/tests/end-to-end-tests/hd_33632Ab.py new file mode 100644 index 00000000..9a193b62 --- /dev/null +++ b/tests/end-to-end-tests/hd_33632Ab.py @@ -0,0 +1,57 @@ +""" +Compare to Hsu+ 2024 (case with HGCA, no RVs, and relative astrometry) +""" + +import os +from orbitize import DATADIR, hipparcos, gaia, read_input, system, priors, sampler + +# the necessary input data for beta Pic is part of the orbitize! example data! +iad_filepath = os.path.join(DATADIR, "H024332.d") +gost_filepath = os.path.join(DATADIR, "gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv") + +# Create the HGCA and helper Hipparcos object +hipparcos_lnprob = hipparcos.HipparcosLogProb(iad_filepath, 24332, 1) +hgca_lnprob = gaia.HGCALogProb(24332, hipparcos_lnprob, gost_filepath) + +# read in relative astrometry +astrometry_filepath = os.path.join(DATADIR, "HD33632Ab.csv") +data_table = read_input.read_file(astrometry_filepath) + +# set up the system, passing in hgca_lnprob and setting it fit dynamical mass +stellar_mass = 1.11 +stellar_mass_err = 0.09 +plx = 37.8953 +plx_err = 0.0263 + +this_system = system.System( + 1, + data_table, + stellar_mass, + plx, + mass_err=stellar_mass_err, + plx_err=plx_err, + fit_secondary_mass=True, + gaia=hgca_lnprob, +) + +# adjust the prior on mass to be uniform between 0 and 0.1 Msol +this_system.sys_priors[this_system.param_idx["m1"]] = priors.LogUniformPrior( + 0, 0.1 +) + +# MCMC parameters +n_temps=20 +n_walkers=1000 +n_threads=20 +total_orbits= n_walkers * 50_000 +burn_steps=10_000 +thin=10 + +total_orbits = 100 * n_walkers + +# create the sampler, run it, and save posteriors +this_sampler = sampler.MCMC(this_system, n_temps, n_walkers, n_threads) + +this_sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + +this_sampler.results.save_results("HD_33632_Ab.hdf5") \ No newline at end of file From dd83a85c67df800eb58354356090a060e8c30395 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 3 Aug 2026 14:17:18 -0700 Subject: [PATCH 06/13] as -> mas --- orbitize/example_data/HD33632Ab.csv | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/orbitize/example_data/HD33632Ab.csv b/orbitize/example_data/HD33632Ab.csv index 3e3fae57..0252e0b3 100644 --- a/orbitize/example_data/HD33632Ab.csv +++ b/orbitize/example_data/HD33632Ab.csv @@ -1,4 +1,4 @@ epoch,object,raoff,raoff_err,decoff,decoff_err,inst -58409.0,1,-0.761,0.005,-0.176,0.004,CHARIS -58423.0,1,-0.753,0.005,-0.178,0.005,NIRC2 -59092.0,1,-0.740,0.005,-0.095,0.003,CHARIS \ No newline at end of file +58409.0,1,-761,5,-176,4,CHARIS +58423.0,1,-753,5,-178,5,NIRC2 +59092.0,1,-740,5,-95,3,CHARIS \ No newline at end of file From 5f19c8579e4faf37fa6bfa502037dfe7ccf4a8d2 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 3 Aug 2026 14:19:27 -0700 Subject: [PATCH 07/13] add periodic save freq & plot --- tests/end-to-end-tests/hd_33632Ab.py | 38 +++++++++++++++++++++------- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/tests/end-to-end-tests/hd_33632Ab.py b/tests/end-to-end-tests/hd_33632Ab.py index 9a193b62..affe88fb 100644 --- a/tests/end-to-end-tests/hd_33632Ab.py +++ b/tests/end-to-end-tests/hd_33632Ab.py @@ -1,9 +1,12 @@ """ Compare to Hsu+ 2024 (case with HGCA, no RVs, and relative astrometry) + +https://iopscience.iop.org/article/10.3847/1538-4357/ad58d3#apjad58d3t5 """ import os -from orbitize import DATADIR, hipparcos, gaia, read_input, system, priors, sampler +from orbitize import DATADIR, hipparcos, gaia, read_input, system, priors, sampler, results +import matplotlib.pyplot as plt # the necessary input data for beta Pic is part of the orbitize! example data! iad_filepath = os.path.join(DATADIR, "H024332.d") @@ -35,9 +38,9 @@ ) # adjust the prior on mass to be uniform between 0 and 0.1 Msol -this_system.sys_priors[this_system.param_idx["m1"]] = priors.LogUniformPrior( - 0, 0.1 -) +# this_system.sys_priors[this_system.param_idx["m1"]] = priors.LogUniformPrior( +# 0, 0.1 +# ) # MCMC parameters n_temps=20 @@ -47,11 +50,28 @@ burn_steps=10_000 thin=10 -total_orbits = 100 * n_walkers +run_fit = True + +if __name__ == '__main__': + + # create the sampler, run it, and save posteriors + this_sampler = sampler.MCMC(this_system, n_temps, n_walkers, n_threads) + + output_filename = "HD_33632_Ab.hdf5" + periodic_save_freq = 5_000 + + if run_fit: + + this_sampler.run_sampler( + total_orbits, burn_steps=burn_steps, thin=thin, periodic_save_freq=periodic_save_freq, + output_filename=output_filename + ) -# create the sampler, run it, and save posteriors -this_sampler = sampler.MCMC(this_system, n_temps, n_walkers, n_threads) + this_sampler.results.save_results(output_filename) -this_sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + myResults = results.Results() + myResults.load_results(output_filename) -this_sampler.results.save_results("HD_33632_Ab.hdf5") \ No newline at end of file + # make corner plot + fig = myResults.plot_corner() + plt.savefig("HD_33632_Ab.png", dpi=250) \ No newline at end of file From 0fba3b78d7a82d280cb8b1daafac7c137d9a8336 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Wed, 5 Aug 2026 15:32:49 -0700 Subject: [PATCH 08/13] add option to load from saved fit --- tests/end-to-end-tests/hd_33632Ab.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/tests/end-to-end-tests/hd_33632Ab.py b/tests/end-to-end-tests/hd_33632Ab.py index affe88fb..fa78e72e 100644 --- a/tests/end-to-end-tests/hd_33632Ab.py +++ b/tests/end-to-end-tests/hd_33632Ab.py @@ -48,16 +48,21 @@ n_threads=20 total_orbits= n_walkers * 50_000 burn_steps=10_000 -thin=10 +thin=100 run_fit = True if __name__ == '__main__': # create the sampler, run it, and save posteriors - this_sampler = sampler.MCMC(this_system, n_temps, n_walkers, n_threads) - output_filename = "HD_33632_Ab.hdf5" + if os.path.exists(output_filename): + prev_result_filename = output_filename + else: + prev_result_filename = None + + this_sampler = sampler.MCMC(this_system, n_temps, n_walkers, n_threads, prev_result_filename=prev_result_filename) + periodic_save_freq = 5_000 if run_fit: From c30d2e8f2d226651c45993934873f7d96643bd0f Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 11 Aug 2026 13:22:57 -0400 Subject: [PATCH 09/13] only create pool if using emcee; ptemcee creates own pool --- orbitize/sampler.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 09c23cc3..c8358b2e 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -993,8 +993,7 @@ def run_sampler( if nsteps <= 0: raise ValueError("Total_orbits must be greater than num_walkers.") - with mp.Pool(processes=self.num_threads) as pool: - if self.use_pt: + if self.use_pt: sampler = ptemcee.Sampler( self.num_walkers, self.num_params, @@ -1006,7 +1005,8 @@ def run_sampler( self.priors, ], ) - else: + else: + with mp.Pool(processes=self.num_threads) as pool: sampler = emcee.EnsembleSampler( self.num_walkers, self.num_params, @@ -1014,7 +1014,7 @@ def run_sampler( pool=pool, kwargs={"include_logp": True}, ) - + print("Starting Burn in") for i, state in enumerate( sampler.sample(self.curr_pos, iterations=burn_steps, thin=thin) From a01a785ddf37b26eea2da0459a350f771b0f9224 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 11 Aug 2026 13:50:26 -0400 Subject: [PATCH 10/13] Revert "only create pool if using emcee; ptemcee creates own pool" This reverts commit c30d2e8f2d226651c45993934873f7d96643bd0f. --- orbitize/sampler.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/orbitize/sampler.py b/orbitize/sampler.py index c8358b2e..09c23cc3 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -993,7 +993,8 @@ def run_sampler( if nsteps <= 0: raise ValueError("Total_orbits must be greater than num_walkers.") - if self.use_pt: + with mp.Pool(processes=self.num_threads) as pool: + if self.use_pt: sampler = ptemcee.Sampler( self.num_walkers, self.num_params, @@ -1005,8 +1006,7 @@ def run_sampler( self.priors, ], ) - else: - with mp.Pool(processes=self.num_threads) as pool: + else: sampler = emcee.EnsembleSampler( self.num_walkers, self.num_params, @@ -1014,7 +1014,7 @@ def run_sampler( pool=pool, kwargs={"include_logp": True}, ) - + print("Starting Burn in") for i, state in enumerate( sampler.sample(self.curr_pos, iterations=burn_steps, thin=thin) From 974702b5a303f9cbd02a335984f5031366704f41 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 11 Aug 2026 13:54:08 -0400 Subject: [PATCH 11/13] pass pool into ptemcee instead of threads --- orbitize/sampler.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 09c23cc3..888f66ab 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -1001,10 +1001,10 @@ def run_sampler( self._logl, orbitize.priors.all_lnpriors, ntemps=self.num_temps, - threads=self.num_threads, logpargs=[ self.priors, ], + pool=pool ) else: sampler = emcee.EnsembleSampler( @@ -1118,7 +1118,7 @@ def run_sampler( self.results.save_results(output_filename) print("Run complete") - # Close pool + if examine_chains: self.examine_chains() From c79e6cd911ce58eb32a739e3ac11cd2dedb690e4 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 1 Sep 2026 12:09:03 -0600 Subject: [PATCH 12/13] Fix quadratic slowdown when using periodic_save_freq in MCMC periodic_save_freq previously re-derived the entire chain history (including an unvectorized per-sample prior recomputation) and rewrote the whole HDF5 output file on every periodic save, making cost grow with total run length instead of just the new chunk. _update_chains_from_sampler now only processes steps since the last save, and Results.save_results appends to resizable HDF5 datasets instead of rewriting the file each time. Co-Authored-By: Claude Sonnet 5 --- orbitize/results.py | 53 +++++++++++++++++++++++++++++--- orbitize/sampler.py | 73 +++++++++++++++++---------------------------- 2 files changed, 77 insertions(+), 49 deletions(-) diff --git a/orbitize/results.py b/orbitize/results.py index fea92e13..37a08ba6 100644 --- a/orbitize/results.py +++ b/orbitize/results.py @@ -48,6 +48,12 @@ def __init__( self.ln_evidence = None self.ln_evidence_err = None + # bookkeeping for save_results(): lets repeated calls for the same + # output file append only the newly added rows instead of rewriting + # the whole file from scratch every time + self._saved_filename = None + self._n_rows_saved = 0 + if self.system is not None: self.tau_ref_epoch = self.system.tau_ref_epoch self.labels = self.system.labels @@ -92,6 +98,24 @@ def add_samples(self, orbital_params, lnlikes, curr_pos=None): if curr_pos is not None: self.curr_pos = curr_pos + def _write_growable_dataset(self, hf, name, data, n_already_saved): + """ + Writes ``data`` to the dataset ``name`` in the open hdf5 file ``hf``. + + If the dataset doesn't exist yet, it's created as resizable so that + future calls can extend it. If it already exists (i.e. a previous + call to ``save_results`` already wrote the first ``n_already_saved`` + rows of ``data`` to this same file), only the new rows beyond that + point are written, instead of rewriting the whole dataset. + """ + if name in hf: + dset = hf[name] + dset.resize(len(data), axis=0) + dset[n_already_saved:] = data[n_already_saved:] + else: + maxshape = (None,) + data.shape[1:] + hf.create_dataset(name, data=data, maxshape=maxshape, chunks=True) + def save_results(self, filename): """ Save results.Results object to an hdf5 file @@ -105,12 +129,26 @@ def save_results(self, filename): ``post``, ``lnlike``, and ``parameter_labels`` are datasets that are members of the root group. + If called repeatedly with the same ``filename`` (e.g. from + ``periodic_save_freq`` during MCMC), only the rows of ``post``/``lnlike`` + that were added since the last call are written to disk, rather than + rewriting the entire accumulated chain every time. + Written: Henry Ngo, 2018 API Update: Sarah Blunt, 2021 """ - hf = h5py.File(filename, 'w') # Creates h5py file object + # if this is a new target file (or the first save), start fresh; + # otherwise reopen the file we've already been writing to and append + if filename != self._saved_filename: + mode = 'w' + self._n_rows_saved = 0 + self._saved_filename = filename + else: + mode = 'a' + + hf = h5py.File(filename, mode) # Creates/opens h5py file object # Add sampler_name as attribute of the root group hf.attrs['sampler_name'] = self.sampler_name @@ -124,15 +162,22 @@ def save_results(self, filename): # Now add post and lnlike from the results object as datasets if self.post is not None: - hf.create_dataset('post', data=self.post) + self._write_growable_dataset(hf, 'post', self.post, self._n_rows_saved) # hf.create_dataset('data', data=self.data) if self.lnlike is not None: - hf.create_dataset('lnlike', data=self.lnlike) + self._write_growable_dataset(hf, 'lnlike', self.lnlike, self._n_rows_saved) + + self._n_rows_saved = len(self.post) if self.post is not None else 0 if self.curr_pos is not None: + if 'curr_pos' in hf: + del hf['curr_pos'] hf.create_dataset("curr_pos", data=self.curr_pos) - self.system.save(hf) + # the system configuration doesn't change between saves, so it only + # needs to be (re)written the first time we save to this file + if mode == 'w': + self.system.save(hf) hf.close() # Closes file object, which writes file to disk diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 888f66ab..17af8fde 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -849,13 +849,19 @@ def _logl(self, params, include_logp=False): return super(MCMC, self)._logl(full_params) + logp - def _update_chains_from_sampler(self, sampler, num_steps=None): + def _update_chains_from_sampler(self, sampler, num_steps=None, start_step=0): """ Updates self.post, self.chain, and self.lnlike from the MCMC sampler Args: sampler (emcee.EnsembleSampler or ptemcee.Sampler): sampler object. num_steps (int): if not None, only stores the first num_steps number of steps + start_step (int): if not 0, only (re)computes steps from this index + onwards, so self.post/self.lnlikes end up covering just the + ``[start_step, num_steps)`` chunk of the chain instead of the + full history. Used during periodic saving so that the (costly, + unvectorized) prior recomputation below doesn't get redone over + steps that were already processed in a previous call. """ if num_steps is None: # use all the steps, grab total number of steps from dimension of chains @@ -866,16 +872,16 @@ def _update_chains_from_sampler(self, sampler, num_steps=None): if self.use_pt: # chain is shape: Ntemp x Nwalkers x Nsteps x Nparams - self.post = sampler.chain[0, :, :num_steps].reshape( + self.post = sampler.chain[0, :, start_step:num_steps].reshape( -1, num_params ) # the reshaping flattens the chain # should also be picking out the lowest temperature logps - self.lnlikes = sampler.loglikelihood[0, :, :num_steps].flatten() - self.lnlikes_alltemps = sampler.loglikelihood[:, :, :num_steps] + self.lnlikes = sampler.loglikelihood[0, :, start_step:num_steps].flatten() + self.lnlikes_alltemps = sampler.loglikelihood[:, :, start_step:num_steps] else: # chain is shape: Nwalkers x Nsteps x Nparams - self.post = sampler.chain[:, :num_steps].reshape(-1, num_params) - self.lnlikes = sampler.lnprobability[:, :num_steps].flatten() + self.post = sampler.chain[:, start_step:num_steps].reshape(-1, num_params) + self.lnlikes = sampler.lnprobability[:, start_step:num_steps].flatten() # convert posterior probability (returned by sampler objects) to likelihood (required by orbitize.results.Results) for i, orb in enumerate(self.post): @@ -1057,62 +1063,39 @@ def run_sampler( if periodic_save_freq is not None: if (i + 1) % periodic_save_freq == 0: # we've completed i+1 steps - self._update_chains_from_sampler(sampler, num_steps=i + 1) - - # figure out what is the new chunk of the chain and corresponding lnlikes that have been computed before last save - # grab the current posterior and lnlikes and reshape them to have the Nwalkers x Nsteps dimension again - post_shape = self.post.shape - curr_chain_shape = ( - self.num_walkers, - post_shape[0] // self.num_walkers, - post_shape[-1], + # only (re)compute the chunk of the chain since the last + # save, not the full history -- avoids redoing the prior + # recomputation over steps that were already saved + self._update_chains_from_sampler( + sampler, num_steps=i + 1, start_step=saved_upto ) - curr_chain = self.post.reshape(curr_chain_shape) - curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2]) - # use the reshaped arrays and find the new steps we computed - curr_chunk = curr_chain[:, saved_upto : i + 1] - curr_chunk = curr_chunk.reshape( - -1, curr_chunk.shape[-1] - ) # flatten nwalkers x nsteps dim - curr_lnlike_chunk = curr_lnlike_chain[ - :, saved_upto : i + 1 - ].flatten() # add this current chunk to the results object (which already has all the previous chunks saved) self.results.add_samples( - curr_chunk, curr_lnlike_chunk, curr_pos=self.curr_pos + self.post, self.lnlikes, curr_pos=self.curr_pos ) self.results.save_results(output_filename) saved_upto = i + 1 print("") - self._update_chains_from_sampler(sampler) if periodic_save_freq is None: - # need to save everything + # nothing has been processed/saved yet; need to do the whole chain + self._update_chains_from_sampler(sampler) self.results.add_samples( self.post, self.lnlikes, curr_pos=self.curr_pos ) elif saved_upto < nsteps: - # just need to save the last few - # same code as above except we just need to grab the last few - post_shape = self.post.shape - curr_chain_shape = ( - self.num_walkers, - post_shape[0] // self.num_walkers, - post_shape[-1], - ) - curr_chain = self.post.reshape(curr_chain_shape) - curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2]) - curr_chunk = curr_chain[:, saved_upto:] - curr_chunk = curr_chunk.reshape( - -1, curr_chunk.shape[-1] - ) # flatten nwalkers x nsteps dim - curr_lnlike_chunk = curr_lnlike_chain[:, saved_upto:].flatten() - + # just need to process and save the leftover chunk since the last save + self._update_chains_from_sampler(sampler, start_step=saved_upto) self.results.add_samples( - curr_chunk, curr_lnlike_chunk, curr_pos=self.curr_pos + self.post, self.lnlikes, curr_pos=self.curr_pos ) + else: + # everything has already been saved; still refresh self.chain, + # self.post, and self.lnlikes to reflect the full run in case + # the caller inspects them after run_sampler() returns + self._update_chains_from_sampler(sampler) if output_filename is not None: self.results.save_results(output_filename) From 58c55f0454f966834e60386214b9b056e1d68ac2 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Fri, 4 Sep 2026 15:12:03 -0700 Subject: [PATCH 13/13] Fix remaining quadratic growth in Results.add_samples during periodic saves Results.add_samples() used np.vstack/np.append, which reallocate and copy the entire accumulated post/lnlike arrays on every call. This meant that even after the earlier fix to periodic_save_freq (only recomputing/writing new chunks), the in-memory accumulation itself was still O(total rows so far) per periodic save, making the MCMC fit progressively slower as it ran. add_samples now writes into a capacity-doubling internal buffer, so appends are amortized O(1) per row instead of O(n) per call. Co-Authored-By: Claude Sonnet 5 --- orbitize/results.py | 63 +++++++++++++++++++++++++++++++++++---------- 1 file changed, 50 insertions(+), 13 deletions(-) diff --git a/orbitize/results.py b/orbitize/results.py index 37a08ba6..b3f5c117 100644 --- a/orbitize/results.py +++ b/orbitize/results.py @@ -54,6 +54,13 @@ def __init__( self._saved_filename = None self._n_rows_saved = 0 + # backing buffers for add_samples(): allow post/lnlike to grow with + # amortized O(1) cost per row instead of reallocating+copying the + # full array on every call (see add_samples() docstring) + self._post_buf = None + self._lnlike_buf = None + self._n_used = 0 + if self.system is not None: self.tau_ref_epoch = self.system.tau_ref_epoch self.labels = self.system.labels @@ -64,36 +71,66 @@ def __init__( self.param_idx = self.system.param_idx self.standard_param_idx = self.system.basis.standard_basis_idx - def add_samples(self, orbital_params, lnlikes, curr_pos=None): + def add_samples(self, orbital_params, lnlikes, curr_pos=None): """ - Add accepted orbits, their likelihoods, and the orbitize version number + Add accepted orbits, their likelihoods, and the orbitize version number to the results Args: - orbital_params (np.array): add sets of orbital params (could be multiple) + orbital_params (np.array): add sets of orbital params (could be multiple) to results lnlike (np.array): add corresponding lnlike values to results - curr_pos (np.array of float): for MCMC only. A multi-D array of the + curr_pos (np.array of float): for MCMC only. A multi-D array of the current walker positions + .. Note:: ``post``/``lnlike`` are backed by an internal buffer that's + over-allocated and grown by doubling, so repeated calls (e.g. from + ``periodic_save_freq`` during MCMC) append in amortized O(1) time + per row instead of reallocating and copying the full accumulated + array on every call. + Written: Henry Ngo, 2018 API Update: Sarah Blunt, 2021 """ - + # Adding the orbitize version number to the results if self.version_number is None: self.version_number = orbitize.__version__ - # If no exisiting results then it is easy - if self.post is None: - self.post = orbital_params - self.lnlike = lnlikes + n_new = len(orbital_params) - # Otherwise, need to append properly - else: - self.post = np.vstack((self.post, orbital_params)) - self.lnlike = np.append(self.lnlike, lnlikes) + # lazily adopt any existing post/lnlike (e.g. set directly via the + # constructor) as the starting buffer + if self._post_buf is None: + if self.post is not None: + self._post_buf = self.post + self._lnlike_buf = self.lnlike + self._n_used = len(self.post) + else: + self._post_buf = np.empty((0, orbital_params.shape[1]), dtype=orbital_params.dtype) + self._lnlike_buf = np.empty(0, dtype=lnlikes.dtype) + self._n_used = 0 + + needed = self._n_used + n_new + capacity = self._post_buf.shape[0] + if needed > capacity: + new_capacity = max(needed, capacity * 2) + + new_post_buf = np.empty((new_capacity,) + self._post_buf.shape[1:], dtype=self._post_buf.dtype) + new_post_buf[: self._n_used] = self._post_buf[: self._n_used] + self._post_buf = new_post_buf + + new_lnlike_buf = np.empty(new_capacity, dtype=self._lnlike_buf.dtype) + new_lnlike_buf[: self._n_used] = self._lnlike_buf[: self._n_used] + self._lnlike_buf = new_lnlike_buf + + self._post_buf[self._n_used : needed] = orbital_params + self._lnlike_buf[self._n_used : needed] = lnlikes + self._n_used = needed + + self.post = self._post_buf[: self._n_used] + self.lnlike = self._lnlike_buf[: self._n_used] if curr_pos is not None: self.curr_pos = curr_pos