From 5a9429f36c1b467dfc19515c05330018748aff74 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Mon, 22 Jul 2024 10:36:25 -0500 Subject: [PATCH 01/44] Passing in arg --- orbitize/kepler.py | 41 ++++++++++++++++++++++++++--------------- orbitize/system.py | 4 ++++ 2 files changed, 30 insertions(+), 15 deletions(-) diff --git a/orbitize/kepler.py b/orbitize/kepler.py index 447bbcc7..9e8c78ea 100644 --- a/orbitize/kepler.py +++ b/orbitize/kepler.py @@ -44,6 +44,29 @@ def tau_to_manom(date, sma, mtot, tau, tau_ref_epoch): return mean_anom +def times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=58849, tolerance=1e-9, max_iter=100, use_c=True, use_gpu=False ): + + n_orbs = np.size(sma) # num sets of input orbital parameters + n_dates = np.size(epochs) # number of dates to compute offsets and vz + + + # Necessary for _calc_ecc_anom, for now + if np.isscalar(epochs): # just in case epochs is given as a scalar + epochs = np.array([epochs]) + ecc_arr = np.tile(ecc, (n_dates, 1)) + + # # compute mean anomaly (size: n_orbs x n_dates) + manom = tau_to_manom(epochs[:, None], sma, mtot, tau, tau_ref_epoch) + # compute eccentric anomalies (size: n_orbs x n_dates) + eanom = _calc_ecc_anom(manom, ecc_arr, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) + + # compute the true anomalies (size: n_orbs x n_dates) + # Note: matrix multiplication makes the shapes work out here and below + tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom)) + + return tanom, eanom + + def calc_orbit( epochs, sma, ecc, inc, aop, pan, tau, plx, mtot, mass_for_Kamp=None, tau_ref_epoch=58849, tolerance=1e-9, @@ -89,26 +112,14 @@ def calc_orbit( Written: Jason Wang, Henry Ngo, 2018 """ - n_orbs = np.size(sma) # num sets of input orbital parameters - n_dates = np.size(epochs) # number of dates to compute offsets and vz - # return planetary RV if `mass_for_Kamp` is not defined + # return planetary RV if `mass_for_Kamp` is not defined if mass_for_Kamp is None: mass_for_Kamp = mtot + ecc - # Necessary for _calc_ecc_anom, for now - if np.isscalar(epochs): # just in case epochs is given as a scalar - epochs = np.array([epochs]) - ecc_arr = np.tile(ecc, (n_dates, 1)) + tanom, eanom = times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=tau_ref_epoch, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) - # # compute mean anomaly (size: n_orbs x n_dates) - manom = tau_to_manom(epochs[:, None], sma, mtot, tau, tau_ref_epoch) - # compute eccentric anomalies (size: n_orbs x n_dates) - eanom = _calc_ecc_anom(manom, ecc_arr, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) - - # compute the true anomalies (size: n_orbs x n_dates) - # Note: matrix multiplication makes the shapes work out here and below - tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom)) # compute 3-D orbital radius of second body (size: n_orbs x n_dates) radius = sma * (1.0 - ecc * np.cos(eanom)) diff --git a/orbitize/system.py b/orbitize/system.py index 1da30d3f..600d8608 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -483,6 +483,10 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): tau_ref_epoch=self.tau_ref_epoch, ) + tanom, eanom = kepler.times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch) + + + # raoff, decoff, vz are scalers if the length of epochs is 1 if len(epochs) == 1: raoff = np.array([raoff]) From d1f285c0245066c26717c223ce26519e03b414d8 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 22 Jul 2024 08:42:30 -0700 Subject: [PATCH 02/44] test commit --- orbitize/kepler.py | 215 ++++++++++++++++++++++++++++++--------------- 1 file changed, 142 insertions(+), 73 deletions(-) diff --git a/orbitize/kepler.py b/orbitize/kepler.py index 9e8c78ea..c3bae2be 100644 --- a/orbitize/kepler.py +++ b/orbitize/kepler.py @@ -1,6 +1,7 @@ """ This module solves for the orbit of the planet given Keplerian parameters. """ + import numpy as np import astropy.units as u import astropy.constants as consts @@ -13,30 +14,29 @@ if cuda_ext: # Configure GPU context for CUDA accelerated compute from orbitize import gpu_context + kep_gpu_ctx = gpu_context.gpu_context() + def tau_to_manom(date, sma, mtot, tau, tau_ref_epoch): """ Gets the mean anomlay - + Args: date (float or np.array): MJD sma (float): semi major axis (AU) mtot (float): total mass (M_sun) tau (float): epoch of periastron, in units of the orbital period tau_ref_epoch (float): reference epoch for tau - + Returns: float or np.array: mean anomaly on that date [0, 2pi) """ - period = np.sqrt( - 4 * np.pi**2.0 * (sma * u.AU)**3 / - (consts.G * (mtot * u.Msun)) - ) + period = np.sqrt(4 * np.pi**2.0 * (sma * u.AU) ** 3 / (consts.G * (mtot * u.Msun))) period = period.to(u.day).value - frac_date = (date - tau_ref_epoch)/period + frac_date = (date - tau_ref_epoch) / period frac_date %= 1 mean_anom = (frac_date - tau) * 2 * np.pi @@ -44,12 +44,25 @@ def tau_to_manom(date, sma, mtot, tau, tau_ref_epoch): return mean_anom -def times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=58849, tolerance=1e-9, max_iter=100, use_c=True, use_gpu=False ): + +def times2trueanom_and_eccanom( + sma, + epochs, + mtot, + ecc, + tau, + tau_ref_epoch=58849, + tolerance=1e-9, + max_iter=100, + use_c=True, + use_gpu=False, +): + + print("hi Farrah!!") n_orbs = np.size(sma) # num sets of input orbital parameters n_dates = np.size(epochs) # number of dates to compute offsets and vz - # Necessary for _calc_ecc_anom, for now if np.isscalar(epochs): # just in case epochs is given as a scalar epochs = np.array([epochs]) @@ -58,21 +71,39 @@ def times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=58849, # # compute mean anomaly (size: n_orbs x n_dates) manom = tau_to_manom(epochs[:, None], sma, mtot, tau, tau_ref_epoch) # compute eccentric anomalies (size: n_orbs x n_dates) - eanom = _calc_ecc_anom(manom, ecc_arr, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) + eanom = _calc_ecc_anom( + manom, + ecc_arr, + tolerance=tolerance, + max_iter=max_iter, + use_c=use_c, + use_gpu=use_gpu, + ) # compute the true anomalies (size: n_orbs x n_dates) # Note: matrix multiplication makes the shapes work out here and below - tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom)) + tanom = 2.0 * np.arctan(np.sqrt((1.0 + ecc) / (1.0 - ecc)) * np.tan(0.5 * eanom)) return tanom, eanom - def calc_orbit( - epochs, sma, ecc, inc, aop, pan, tau, plx, mtot, mass_for_Kamp=None, tau_ref_epoch=58849, tolerance=1e-9, - max_iter=100, use_c=True, use_gpu=False + epochs, + sma, + ecc, + inc, + aop, + pan, + tau, + plx, + mtot, + mass_for_Kamp=None, + tau_ref_epoch=58849, + tolerance=1e-9, + max_iter=100, + use_c=True, + use_gpu=False, ): - """ Returns the separation and radial velocity of the body given array of orbital parameters (size n_orbs) at given epochs (array of size n_dates) @@ -113,20 +144,31 @@ def calc_orbit( Written: Jason Wang, Henry Ngo, 2018 """ - # return planetary RV if `mass_for_Kamp` is not defined + # return planetary RV if `mass_for_Kamp` is not defined if mass_for_Kamp is None: mass_for_Kamp = mtot ecc - tanom, eanom = times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=tau_ref_epoch, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) + tanom, eanom = times2trueanom_and_eccanom( + sma, + epochs, + mtot, + ecc, + tau, + tau_ref_epoch=tau_ref_epoch, + tolerance=tolerance, + max_iter=max_iter, + use_c=use_c, + use_gpu=use_gpu, + ) # compute 3-D orbital radius of second body (size: n_orbs x n_dates) radius = sma * (1.0 - ecc * np.cos(eanom)) # compute ra/dec offsets (size: n_orbs x n_dates) # math from James Graham. Lots of trig - c2i2 = np.cos(0.5*inc)**2 - s2i2 = np.sin(0.5*inc)**2 + c2i2 = np.cos(0.5 * inc) ** 2 + s2i2 = np.sin(0.5 * inc) ** 2 arg1 = tanom + aop + pan arg2 = tanom + aop - pan c1 = np.cos(arg1) @@ -135,39 +177,46 @@ def calc_orbit( s2 = np.sin(arg2) # updated sign convention for Green Eq. 19.4-19.7 - raoff = radius * (c2i2*s1 - s2i2*s2) * plx - deoff = radius * (c2i2*c1 + s2i2*c2) * plx + raoff = radius * (c2i2 * s1 - s2i2 * s2) * plx + deoff = radius * (c2i2 * c1 + s2i2 * c2) * plx # compute the radial velocity (vz) of the body (size: n_orbs x n_dates) # first comptue the RV semi-amplitude (size: n_orbs x n_dates) - Kv = np.sqrt(consts.G / (1.0 - ecc**2)) * (mass_for_Kamp * u.Msun * - np.sin(inc)) / np.sqrt(mtot * u.Msun) / np.sqrt(sma * u.au) + Kv = ( + np.sqrt(consts.G / (1.0 - ecc**2)) + * (mass_for_Kamp * u.Msun * np.sin(inc)) + / np.sqrt(mtot * u.Msun) + / np.sqrt(sma * u.au) + ) # Convert to km/s - Kv = Kv.to(u.km/u.s) + Kv = Kv.to(u.km / u.s) # compute the vz - vz = Kv.value * (ecc*np.cos(aop) + np.cos(aop + tanom)) + vz = Kv.value * (ecc * np.cos(aop) + np.cos(aop + tanom)) # Squeeze out extra dimension (useful if n_orbs = 1, does nothing if n_orbs > 1) vz = np.squeeze(vz)[()] return raoff, deoff, vz -def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gpu=False): + +def _calc_ecc_anom( + manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gpu=False +): """ - Computes the eccentric anomaly from the mean anomlay. - Code from Rob De Rosa's orbit solver (e < 0.95 use Newton, e >= 0.95 use Mikkola) + Computes the eccentric anomaly from the mean anomlay. + Code from Rob De Rosa's orbit solver (e < 0.95 use Newton, e >= 0.95 use Mikkola) - Args: - manom (float/np.array): mean anomaly, either a scalar or np.array of any shape - ecc (float/np.array): eccentricity, either a scalar or np.array of the same shape as manom - tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. - max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. - use_c (bool, optional): Use the C solver if configured. Defaults to False - use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False + Args: + manom (float/np.array): mean anomaly, either a scalar or np.array of any shape + ecc (float/np.array): eccentricity, either a scalar or np.array of the same shape as manom + tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. + max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. + use_c (bool, optional): Use the C solver if configured. Defaults to False + use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False -Return: - eanom (float/np.array): eccentric anomalies, same shape as manom + Return: + eanom (float/np.array): eccentric anomalies, same shape as manom - Written: Jason Wang, 2018 + Written: Jason Wang, 2018 """ if np.isscalar(ecc) or (np.shape(manom) == np.shape(ecc)): @@ -177,7 +226,7 @@ def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gp # If manom is a scalar, make it into a one-element array if np.isscalar(manom): - manom = np.array((manom, )) + manom = np.array((manom,)) # If ecc is a scalar, make it the same shape as manom if np.isscalar(ecc): @@ -197,19 +246,26 @@ def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gp # Now low eccentricities ind_low = np.where(~ecc_zero & ecc_low) - if len(ind_low[0]) > 0: - eanom[ind_low] = _newton_solver_wrapper(manom[ind_low], ecc[ind_low], tolerance, max_iter, use_c, use_gpu) - + if len(ind_low[0]) > 0: + eanom[ind_low] = _newton_solver_wrapper( + manom[ind_low], ecc[ind_low], tolerance, max_iter, use_c, use_gpu + ) + # Now high eccentricities - ind_high = np.where(~ecc_zero & ~ecc_low | (eanom == -1)) # The C and CUDA solvers return the unphysical value -1 if they fail to converge - if len(ind_high[0]) > 0: - eanom[ind_high] = _mikkola_solver_wrapper(manom[ind_high], ecc[ind_high], use_c, use_gpu) + ind_high = np.where( + ~ecc_zero & ~ecc_low | (eanom == -1) + ) # The C and CUDA solvers return the unphysical value -1 if they fail to converge + if len(ind_high[0]) > 0: + eanom[ind_high] = _mikkola_solver_wrapper( + manom[ind_high], ecc[ind_high], use_c, use_gpu + ) return np.squeeze(eanom)[()] + def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False, use_gpu=False): """ - Wrapper for the various (Python, C, CUDA) implementations of the Newton-Raphson solver + Wrapper for the various (Python, C, CUDA) implementations of the Newton-Raphson solver for eccentric anomaly. Args: @@ -224,18 +280,21 @@ def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False, use_gpu Written: Devin Cody, 2021 """ eanom = np.empty_like(manom) - + if cuda_ext and use_gpu: # the CUDA solver returns eanom = -1 if it doesnt converge after max_iter iterations eanom = _CUDA_newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) elif cext and use_c: # the C solver returns eanom = -1 if it doesnt converge after max_iter iterations - eanom = _kepler._c_newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) + eanom = _kepler._c_newton_solver( + manom, ecc, tolerance=tolerance, max_iter=max_iter + ) else: eanom = _newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) return eanom + def _newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): """ Newton-Raphson solver for eccentric anomaly. @@ -243,11 +302,11 @@ def _newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): Args: manom (np.array): array of mean anomalies ecc (np.array): array of eccentricities - tolerance (float, optional): absolute tolerance of iterative computation. + tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. - max_iter (int, optional): maximum number of iterations before switching. + max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. - eanom0 (np.array): array of first guess for eccentric anomaly, same + eanom0 (np.array): array of first guess for eccentric anomaly, same shape as manom (optional) @@ -273,23 +332,29 @@ def _newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): abs_diff = np.abs(diff) ind = np.where(abs_diff > tolerance) niter = 0 - while ((ind[0].size > 0) and (niter <= max_iter)): + while (ind[0].size > 0) and (niter <= max_iter): eanom[ind] -= diff[ind] # If it hasn't converged after half the iterations are done, try starting from pi - if niter == (max_iter//2): + if niter == (max_iter // 2): eanom[ind] = np.pi - diff[ind] = (eanom[ind] - (ecc[ind] * np.sin(eanom[ind])) - manom[ind]) / \ - (1.0 - (ecc[ind] * np.cos(eanom[ind]))) + diff[ind] = (eanom[ind] - (ecc[ind] * np.sin(eanom[ind])) - manom[ind]) / ( + 1.0 - (ecc[ind] * np.cos(eanom[ind])) + ) abs_diff[ind] = np.abs(diff[ind]) ind = np.where(abs_diff > tolerance) niter += 1 if niter >= max_iter: - print(manom[ind], eanom[ind], diff[ind], ecc[ind], '> {} iter.'.format(max_iter)) - eanom[ind] = _mikkola_solver_wrapper(manom[ind], ecc[ind]) # Send remaining orbits to the analytical version, this has not happened yet... + print( + manom[ind], eanom[ind], diff[ind], ecc[ind], "> {} iter.".format(max_iter) + ) + eanom[ind] = _mikkola_solver_wrapper( + manom[ind], ecc[ind] + ) # Send remaining orbits to the analytical version, this has not happened yet... return eanom + def _CUDA_newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): """ Helper function for calling the CUDA implementation of the Newton-Raphson solver for eccentric anomaly. @@ -309,16 +374,17 @@ def _CUDA_newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): manom = np.asarray(manom) ecc = np.asarray(ecc) eanom = np.empty_like(manom) - tolerance = np.asarray(tolerance, dtype = np.float64) + tolerance = np.asarray(tolerance, dtype=np.float64) max_iter = np.asarray(max_iter) - + kep_gpu_ctx.newton(manom, ecc, eanom, eanom0, tolerance, max_iter) return eanom + def _mikkola_solver_wrapper(manom, ecc, use_c=False, use_gpu=False): """ - Wrapper for the various (Python, C, CUDA) implementations of Analtyical Mikkola solver + Wrapper for the various (Python, C, CUDA) implementations of Analtyical Mikkola solver Args: manom (np.array): array of mean anomalies between 0 and 2pi @@ -364,26 +430,29 @@ def _mikkola_solver(manom, ecc): beta = (0.5 * manom) / ((4.0 * ecc) + 0.5) aux = np.sqrt(beta**2.0 + alpha**3.0) - z = np.abs(beta + aux)**(1.0/3.0) + z = np.abs(beta + aux) ** (1.0 / 3.0) - s0 = z - (alpha/z) - s1 = s0 - (0.078*(s0**5.0)) / (1.0 + ecc) - e0 = manom + (ecc * (3.0*s1 - 4.0*(s1**3.0))) + s0 = z - (alpha / z) + s1 = s0 - (0.078 * (s0**5.0)) / (1.0 + ecc) + e0 = manom + (ecc * (3.0 * s1 - 4.0 * (s1**3.0))) se0 = np.sin(e0) ce0 = np.cos(e0) - f = e0-ecc*se0-manom - f1 = 1.0-ecc*ce0 - f2 = ecc*se0 - f3 = ecc*ce0 + f = e0 - ecc * se0 - manom + f1 = 1.0 - ecc * ce0 + f2 = ecc * se0 + f3 = ecc * ce0 f4 = -f2 - u1 = -f/f1 - u2 = -f/(f1+0.5*f2*u1) - u3 = -f/(f1+0.5*f2*u2+(1.0/6.0)*f3*u2*u2) - u4 = -f/(f1+0.5*f2*u3+(1.0/6.0)*f3*u3*u3+(1.0/24.0)*f4*(u3**3.0)) + u1 = -f / f1 + u2 = -f / (f1 + 0.5 * f2 * u1) + u3 = -f / (f1 + 0.5 * f2 * u2 + (1.0 / 6.0) * f3 * u2 * u2) + u4 = -f / ( + f1 + 0.5 * f2 * u3 + (1.0 / 6.0) * f3 * u3 * u3 + (1.0 / 24.0) * f4 * (u3**3.0) + ) + + return e0 + u4 - return (e0 + u4) def _CUDA_mikkola_solver(manom, ecc): """ From a0b603ffb412ef0d7a2e9ef5657d21ddebfd6b94 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 22 Jul 2024 08:58:45 -0700 Subject: [PATCH 03/44] add minimal unit test for reflected-light functionality --- tests/test_brightness.py | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 tests/test_brightness.py diff --git a/tests/test_brightness.py b/tests/test_brightness.py new file mode 100644 index 00000000..f545ec3e --- /dev/null +++ b/tests/test_brightness.py @@ -0,0 +1,40 @@ +""" +Tests for the reflected-light calculation in system.compute_all_orbits +""" + +from orbitize import system +from orbitize import DATADIR, read_input +import os +import numpy as np + + +def test_brightness_calculation(): + + num_secondary_bodies = 1 + input_file = os.path.join(DATADIR, "GJ504.csv") + data_table = read_input.read_file(input_file) + system_mass = 1.47 + plx = 24.30 + + test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) + + params = np.array( + [ + 7.2774010e01, + 4.1116819e-02, + 5.6322372e-01, + 3.5251172e00, + 4.2904768e00, + 9.4234377e-02, + 4.5418411e01, + 1.4317369e-03, + 5.6322372e-01, + 3.1016846e00, + 4.2904768e00, + 3.4033456e-01, + 2.4589758e01, + 1.4849439e00, + ] + ) + + ra, dec, vz, brightness = test_system.compute_all_orbits(params) From 04e86384057bbf2bb5b19ba69e26a842866aca59 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 22 Jul 2024 08:59:31 -0700 Subject: [PATCH 04/44] add if name == main to test_brightness --- tests/test_brightness.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index f545ec3e..e3fdfedb 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -38,3 +38,7 @@ def test_brightness_calculation(): ) ra, dec, vz, brightness = test_system.compute_all_orbits(params) + + +if __name__ == "__main__": + test_brightness_calculation() From 2d9129b4c88df1561d203718c40a42ff7b217095 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Tue, 23 Jul 2024 14:26:02 -0500 Subject: [PATCH 05/44] Called times2trueanom_and_eccanom --- orbitize/system.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orbitize/system.py b/orbitize/system.py index 600d8608..bc7ae6aa 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -483,7 +483,7 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): tau_ref_epoch=self.tau_ref_epoch, ) - tanom, eanom = kepler.times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch) + tanom, eanom = kepler.times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=58849, tolerance=1e-9, max_iter=100, use_c=True, use_gpu=False,) From 017c8ce30b1d35410cf0cf74a8454417651796c2 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Wed, 24 Jul 2024 09:57:51 -0700 Subject: [PATCH 06/44] modify read_input to read in brightness vals --- orbitize/read_input.py | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/orbitize/read_input.py b/orbitize/read_input.py index f37f5e41..69e2f0fe 100644 --- a/orbitize/read_input.py +++ b/orbitize/read_input.py @@ -103,8 +103,9 @@ def read_file(filename): "quant12_corr", "quant_type", "instrument", + "brightness", ), - dtype=(float, int, float, float, float, float, float, "S5", "S5"), + dtype=(float, int, float, float, float, float, float, "S5", "S5", float), ) # read file @@ -179,6 +180,10 @@ def read_file(filename): have_seppacorr = np.zeros( num_measurements, dtype=bool ) # zeros are False + if "brightness" in input_table.columns: + have_brightness = ~input_table["brightness"].mask + else: + have_brightness = np.zeros(num_measurements, dtype=bool) if "rv" in input_table.columns: have_rv = ~input_table["rv"].mask else: @@ -231,11 +236,14 @@ def read_file(filename): else: have_rv = np.zeros(num_measurements, dtype=bool) # zeros are False - # Rob: not sure if we need this but adding just in case if "instrument" in input_table.columns: have_inst = np.ones(num_measurements, dtype=bool) else: have_inst = np.zeros(num_measurements, dtype=bool) + if "brightness" in input_table.columns: + have_brightness = np.ones(num_measurements, dtype=bool) + else: + have_brightness = np.zeros(num_measurements, dtype=bool) # orbitize! backwards compatability since we added new columns, some old data formats may not have them # fill in with default values @@ -243,6 +251,12 @@ def read_file(filename): if "quant12_corr" not in input_table.keys(): default_corrs = np.nan * np.ones(len(input_table)) input_table.add_column(default_corrs, name="quant12_corr") + if "brightness" not in input_table.keys(): + default_brightness = np.nan * np.ones(len(input_table)) + input_table.add_column(default_brightness, name="brightness") + have_brightness = np.zeros(num_measurements, dtype=bool) + else: + have_brightness = np.ones(num_measurements, dtype=bool) if "instrument" not in input_table.keys(): default_insts = [] for this_quant_type in input_table["quant_type"]: @@ -275,6 +289,10 @@ def read_file(filename): MJD = row["epoch"] - 2400000.5 else: MJD = row["epoch"] + if have_brightness[index]: + brightness = row["brightness"] + else: + brightness = None # check that "object" is an integer (instead of ABC/bcd) if not isinstance(row["object"], (int, np.int32, np.int64)): @@ -294,6 +312,7 @@ def read_file(filename): None, row["quant_type"], row["instrument"], + None, ] ) @@ -322,6 +341,7 @@ def read_file(filename): quant12_corr, row["quant_type"], row["instrument"], + brightness, ] ) else: # catch wrong formats @@ -332,6 +352,7 @@ def read_file(filename): ) else: # When not in orbitize style + if have_ra[index] and have_dec[index]: # check if there's a covariance term if have_radeccorr[index]: @@ -362,6 +383,7 @@ def read_file(filename): quant12_corr, "radec", this_inst, + brightness, ] ) @@ -395,6 +417,7 @@ def read_file(filename): quant12_corr, "seppa", this_inst, + brightness, ] ) @@ -411,6 +434,7 @@ def read_file(filename): None, "rv", row["instrument"], + brightness, ] ) else: @@ -426,6 +450,7 @@ def read_file(filename): None, "rv", "defrv", + brightness, ] ) From 03ea9844904c291d4cda2f05e58dc6d1f2b3109c Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Wed, 24 Jul 2024 09:58:45 -0700 Subject: [PATCH 07/44] update todos in test_brightness --- tests/test_brightness.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index e3fdfedb..a56daf95 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -11,8 +11,14 @@ def test_brightness_calculation(): num_secondary_bodies = 1 + + # TODO (sarah): change dataset to one where we can see brightness variations input_file = os.path.join(DATADIR, "GJ504.csv") + # input_file = os.path.join(DATADIR, "new_dataset.csv") data_table = read_input.read_file(input_file) + + times = data_table["epoch"].value + system_mass = 1.47 plx = 24.30 @@ -37,8 +43,18 @@ def test_brightness_calculation(): ] ) + print(test_system.param_idx) + ra, dec, vz, brightness = test_system.compute_all_orbits(params) + # TODO (farrah): make plot of brightness vs time + + +def test_read_input_with_brightness(): + + # TODO (farrah): use code above as inspiration to read in a csv file with a brightness column + if __name__ == "__main__": test_brightness_calculation() + test_read_input_with_brightness() From 64045dd54b3830b11e5e40c8a828f614075291d2 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Wed, 24 Jul 2024 10:09:36 -0700 Subject: [PATCH 08/44] test system for brightness calc: gj 504 -> beta pic --- tests/test_brightness.py | 33 ++++++++++++++------------------- 1 file changed, 14 insertions(+), 19 deletions(-) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index a56daf95..1a7afd52 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -12,9 +12,8 @@ def test_brightness_calculation(): num_secondary_bodies = 1 - # TODO (sarah): change dataset to one where we can see brightness variations - input_file = os.path.join(DATADIR, "GJ504.csv") - # input_file = os.path.join(DATADIR, "new_dataset.csv") + # input_file = os.path.join(DATADIR, "GJ504.csv") + input_file = os.path.join(DATADIR, "betaPic.csv") data_table = read_input.read_file(input_file) times = data_table["epoch"].value @@ -24,27 +23,21 @@ def test_brightness_calculation(): test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) + print(test_system.param_idx) + params = np.array( [ - 7.2774010e01, - 4.1116819e-02, - 5.6322372e-01, - 3.5251172e00, - 4.2904768e00, - 9.4234377e-02, - 4.5418411e01, - 1.4317369e-03, - 5.6322372e-01, - 3.1016846e00, - 4.2904768e00, - 3.4033456e-01, - 2.4589758e01, - 1.4849439e00, + 10.0, + 0.1, + np.radians(89), + np.radians(21), + np.radians(31), + 0.0, # note: I didn't convert tau here, just picked random number + 51.5, + 1.75, ] ) - print(test_system.param_idx) - ra, dec, vz, brightness = test_system.compute_all_orbits(params) # TODO (farrah): make plot of brightness vs time @@ -54,6 +47,8 @@ def test_read_input_with_brightness(): # TODO (farrah): use code above as inspiration to read in a csv file with a brightness column + print("hello! :D ") + if __name__ == "__main__": test_brightness_calculation() From 9b9055e82120ce3b4795e102682aa138b75c51cf Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Mon, 29 Jul 2024 16:57:20 -0500 Subject: [PATCH 09/44] Added Brightness_calculation --- orbitize/system.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/orbitize/system.py b/orbitize/system.py index bc7ae6aa..7dce9d95 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -485,6 +485,22 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): tanom, eanom = kepler.times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=58849, tolerance=1e-9, max_iter=100, use_c=True, use_gpu=False,) + def brightness_calculation(sma, tanom, ecc = 0.75, inc = np.radians(30), aop = np.radians(120), plx = np.radians(65)): + + R = (sma*(1-ecc**2))/(1+ecc*np.cos(tanom)) + + z = (R)*(-np.cos(aop)*np.sin(inc)*np.sin(tanom)-np.cos(tanom)*np.sin(inc)*np.sin(aop)) + + #ro = (z**2+((raoff**2+decoff**2)/plx**2))**0.5 + + B = math.arctan2(-R, z)+ math.pi + + Alpha = (1/math.pi)*(np.sin(B)+(math.pi-B)*np.cos(B)) + + Albedo = 0.5 + brightness = Albedo*Alpha/R**2 + + return brightness # raoff, decoff, vz are scalers if the length of epochs is 1 From 81f8388ab6ff7dcace73dd5e442b5d1c9273b334 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 30 Jul 2024 11:33:30 -0700 Subject: [PATCH 10/44] fixed bug returning brightness --- orbitize/system.py | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index 7dce9d95..9e7aef10 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -483,25 +483,20 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): tau_ref_epoch=self.tau_ref_epoch, ) - tanom, eanom = kepler.times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=58849, tolerance=1e-9, max_iter=100, use_c=True, use_gpu=False,) + tanom, eanom = kepler.times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=self.tau_ref_epoch) - def brightness_calculation(sma, tanom, ecc = 0.75, inc = np.radians(30), aop = np.radians(120), plx = np.radians(65)): - R = (sma*(1-ecc**2))/(1+ecc*np.cos(tanom)) + R = (sma*(1-ecc**2))/(1+ecc*np.cos(tanom)) - z = (R)*(-np.cos(aop)*np.sin(inc)*np.sin(tanom)-np.cos(tanom)*np.sin(inc)*np.sin(aop)) - - #ro = (z**2+((raoff**2+decoff**2)/plx**2))**0.5 - - B = math.arctan2(-R, z)+ math.pi + z = (R)*(-np.cos(argp)*np.sin(inc)*np.sin(tanom)-np.cos(tanom)*np.sin(inc)*np.sin(argp)) + + B = np.arctan2(-R, z)+ np.pi - Alpha = (1/math.pi)*(np.sin(B)+(math.pi-B)*np.cos(B)) + Alpha = (1/np.pi)*(np.sin(B)+(np.pi-B)*np.cos(B)) - Albedo = 0.5 - brightness = Albedo*Alpha/R**2 + Albedo = 0.5 + brightness = Albedo*Alpha/R**2 - return brightness - # raoff, decoff, vz are scalers if the length of epochs is 1 if len(epochs) == 1: @@ -584,7 +579,7 @@ def brightness_calculation(sma, tanom, ecc = 0.75, inc = np.radians(30), aop = n else: return raoff, deoff, vz else: - return raoff, deoff, vz + return raoff, deoff, vz, brightness def compute_model(self, params_arr, use_rebound=False): """ From ec0752060f4fbe718dcba5714591da2cc50dd392 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Mon, 5 Aug 2024 10:21:12 -0500 Subject: [PATCH 11/44] Plot of betaPicB data --- tests/test_brightness.py | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 1a7afd52..692d5279 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -1,18 +1,17 @@ """ Tests for the reflected-light calculation in system.compute_all_orbits """ - from orbitize import system from orbitize import DATADIR, read_input import os import numpy as np - +import matplotlib.pyplot as plt def test_brightness_calculation(): num_secondary_bodies = 1 - # input_file = os.path.join(DATADIR, "GJ504.csv") + #input_file = os.path.join(DATADIR, "GJ504.csv") input_file = os.path.join(DATADIR, "betaPic.csv") data_table = read_input.read_file(input_file) @@ -39,16 +38,33 @@ def test_brightness_calculation(): ) ra, dec, vz, brightness = test_system.compute_all_orbits(params) - # TODO (farrah): make plot of brightness vs time + plt.figure() + plt.scatter(times, brightness) + plt.xlabel("Time [dy]") + plt.ylabel("Brightness") + plt.savefig("Test_brightness.png") + def test_read_input_with_brightness(): # TODO (farrah): use code above as inspiration to read in a csv file with a brightness column + num_secondary_bodies = 1 - print("hello! :D ") + # input_file = os.path.join(DATADIR, "GJ504.csv") + input_file = os.path.join(DATADIR, "betaPic.csv") + + data_table = read_input.read_file(input_file) + times = data_table["epoch"].value + brightness_values = data_table["brightness"].value + + +# Do we need the rest of this? since the values for time and brightness are given + + print("hello! :D ") + if __name__ == "__main__": test_brightness_calculation() From f308ed4f518c97e600d7adea9c686a97403db5fe Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Mon, 5 Aug 2024 10:22:33 -0500 Subject: [PATCH 12/44] test --- tests/test_brightness.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 692d5279..37288a76 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -69,3 +69,5 @@ def test_read_input_with_brightness(): if __name__ == "__main__": test_brightness_calculation() test_read_input_with_brightness() + + #Test commit From 4c3f6cf10d1613d6d0595804a93baf62616dfba2 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 5 Aug 2024 09:03:23 -0700 Subject: [PATCH 13/44] add brightness to return stms of comp_all_orbits & todo for farrah --- orbitize/system.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/orbitize/system.py b/orbitize/system.py index 9e7aef10..14a92b66 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -614,7 +614,7 @@ def compute_model(self, params_arr, use_rebound=False): standard_params_arr, comp_rebound=True ) else: - raoff, decoff, vz = self.compute_all_orbits(standard_params_arr) + raoff, decoff, vz, brightness = self.compute_all_orbits(standard_params_arr) if len(standard_params_arr.shape) == 1: n_orbits = 1 @@ -666,6 +666,10 @@ def compute_model(self, params_arr, use_rebound=False): model[self.rv[body_num], 0] = vz[self.rv[body_num], body_num, :] model[self.rv[body_num], 1] = np.nan + # TODO (farrah): add brightness to model here (use RV block above as template) + # (assume self.brightness is array of indices of epochs with brightness measurements) + + # if we have abs astrometry measurements in the input file (i.e. not # from Hipparcos or Gaia), add the parallactic & proper motion here by # calling AbsAstrom compute_model From 37bbb193a3e689534d5e059ea7d881c13f84a21b Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 5 Aug 2024 09:03:41 -0700 Subject: [PATCH 14/44] add test_compute_posteriors --- tests/test_brightness.py | 58 +++++++++++++++++++++++++++++++++------- 1 file changed, 48 insertions(+), 10 deletions(-) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 37288a76..9c5ca348 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -1,17 +1,19 @@ """ Tests for the reflected-light calculation in system.compute_all_orbits """ -from orbitize import system + +from orbitize import system, sampler from orbitize import DATADIR, read_input import os import numpy as np import matplotlib.pyplot as plt + def test_brightness_calculation(): num_secondary_bodies = 1 - #input_file = os.path.join(DATADIR, "GJ504.csv") + # input_file = os.path.join(DATADIR, "GJ504.csv") input_file = os.path.join(DATADIR, "betaPic.csv") data_table = read_input.read_file(input_file) @@ -40,13 +42,13 @@ def test_brightness_calculation(): ra, dec, vz, brightness = test_system.compute_all_orbits(params) # TODO (farrah): make plot of brightness vs time - plt.figure() plt.scatter(times, brightness) plt.xlabel("Time [dy]") plt.ylabel("Brightness") plt.savefig("Test_brightness.png") + def test_read_input_with_brightness(): # TODO (farrah): use code above as inspiration to read in a csv file with a brightness column @@ -54,20 +56,56 @@ def test_read_input_with_brightness(): # input_file = os.path.join(DATADIR, "GJ504.csv") input_file = os.path.join(DATADIR, "betaPic.csv") - + data_table = read_input.read_file(input_file) times = data_table["epoch"].value brightness_values = data_table["brightness"].value + # Do we need the rest of this? since the values for time and brightness are given -# Do we need the rest of this? since the values for time and brightness are given - print("hello! :D ") - + + +def test_compute_posteriors(): + + num_secondary_bodies = 1 + + input_file = os.path.join(DATADIR, "GJ504.csv") + data_table = read_input.read_file(input_file) + + system_mass = 1.47 + plx = 24.30 + + test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) + + params_arr = np.array( + [ + 10.0, + 0.1, + np.radians(89), + np.radians(21), + np.radians(31), + 0.0, # note: I didn't convert tau here, just picked random number + 51.5, + 1.75, + ] + ) + + all_orbits = test_system.compute_all_orbits(params_arr) + print(all_orbits) + + model = test_system.compute_model(params_arr) + # print(model) + + # test_mcmc = sampler.MCMC(test_system, 1, 50, num_threads=1) + + # test_mcmc.run_sampler(10) + if __name__ == "__main__": - test_brightness_calculation() - test_read_input_with_brightness() + # test_brightness_calculation() + # test_read_input_with_brightness() + test_compute_posteriors() - #Test commit + # Test commit From 9dcd884eb62daad5424f0eb3f2c5e6361e3e695b Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Mon, 5 Aug 2024 11:30:30 -0500 Subject: [PATCH 15/44] Added brightness as a model prediction --- orbitize/system.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/orbitize/system.py b/orbitize/system.py index 14a92b66..c016b4ba 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -667,6 +667,10 @@ def compute_model(self, params_arr, use_rebound=False): model[self.rv[body_num], 1] = np.nan # TODO (farrah): add brightness to model here (use RV block above as template) + # Brightness + if len(self.brightness[body_num]) > 0: + model[self.brightness[body_num], 0] = brightness[self.brightness[body_num], body_num, :] + model[self.brightness[body_num], 1] = np.nan # (assume self.brightness is array of indices of epochs with brightness measurements) From 222dec72983262805a98d9e26a334804f4618570 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 19 Aug 2024 10:50:57 -0700 Subject: [PATCH 16/44] add pretty visual to test for farrah --- tests/test_brightness.py | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 9c5ca348..464dec53 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -81,21 +81,29 @@ def test_compute_posteriors(): params_arr = np.array( [ - 10.0, - 0.1, - np.radians(89), - np.radians(21), - np.radians(31), - 0.0, # note: I didn't convert tau here, just picked random number - 51.5, - 1.75, + 10.0, # sma + 0.3, # ecc + np.radians(0), # inc + np.radians(45), # aop + np.radians(90), # pan + 0.0, # tau + 51.5, # plx + 1.75, # stellar maxx ] ) + epochs = np.linspace(0, 365*30, int(1e3)) + ra, dec, vz, brightness = test_system.compute_all_orbits(params_arr, epochs=epochs) + + fig, ax = plt.subplots(2, 1, figsize=(5,10)) + + + ax[0].scatter(epochs, brightness, color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) + ax[1].scatter(ra[:,1,:], dec[:,1,:], color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) - all_orbits = test_system.compute_all_orbits(params_arr) - print(all_orbits) + ax[1].axis('equal') + plt.savefig('visual4farrah.png') - model = test_system.compute_model(params_arr) + # model = test_system.compute_model(params_arr) # print(model) # test_mcmc = sampler.MCMC(test_system, 1, 50, num_threads=1) From 1863981c0043a5e937e7982acb9c54c75ac0c3a8 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 17 Sep 2024 16:11:47 -0700 Subject: [PATCH 17/44] add todo for farrah in read input test --- .../example_data/reflected_light_example.csv | 10 ++++++++++ orbitize/read_input.py | 2 +- tests/test_brightness.py | 17 +++++++---------- 3 files changed, 18 insertions(+), 11 deletions(-) create mode 100644 orbitize/example_data/reflected_light_example.csv diff --git a/orbitize/example_data/reflected_light_example.csv b/orbitize/example_data/reflected_light_example.csv new file mode 100644 index 00000000..5a17cc4f --- /dev/null +++ b/orbitize/example_data/reflected_light_example.csv @@ -0,0 +1,10 @@ +epoch,object,raoff,decoff,raoff_err,decoff_err,brightness +57298,1,253.72,92.35,2.98,2.85,0.5 +57606,1,236.63,127.94,9.77,9.18,0.5 +57645,1,234.52,123.39,1.79,1.03,0.5 +57946,1,210.76,152.09,1.94,1.88,0.5 +58276,1,167.49,180.87,1.61,16.97,0.5 +58287,1,177.67,174.6,1.67,1.67,0.5 +58365,1,165.7,185.33,3.28,3.66,0.5 +58368,1,170.38,185.94,2.52,2.74,0.5 +58414,1,161.64,176.21,13.6,14.31, \ No newline at end of file diff --git a/orbitize/read_input.py b/orbitize/read_input.py index 69e2f0fe..9927c706 100644 --- a/orbitize/read_input.py +++ b/orbitize/read_input.py @@ -292,7 +292,7 @@ def read_file(filename): if have_brightness[index]: brightness = row["brightness"] else: - brightness = None + brightness = np.nan # check that "object" is an integer (instead of ABC/bcd) if not isinstance(row["object"], (int, np.int32, np.int64)): diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 464dec53..48dd1e1f 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -40,7 +40,6 @@ def test_brightness_calculation(): ) ra, dec, vz, brightness = test_system.compute_all_orbits(params) - # TODO (farrah): make plot of brightness vs time plt.figure() plt.scatter(times, brightness) @@ -51,20 +50,20 @@ def test_brightness_calculation(): def test_read_input_with_brightness(): - # TODO (farrah): use code above as inspiration to read in a csv file with a brightness column num_secondary_bodies = 1 - # input_file = os.path.join(DATADIR, "GJ504.csv") - input_file = os.path.join(DATADIR, "betaPic.csv") + input_file = os.path.join(DATADIR, "reflected_light_example.csv") data_table = read_input.read_file(input_file) times = data_table["epoch"].value brightness_values = data_table["brightness"].value - # Do we need the rest of this? since the values for time and brightness are given + print(data_table) - print("hello! :D ") + # TODO (Farrah): add a test that asserts the brightness column of the data table is + # what you expect (hint: check in the reflected_light_example.csv to see what + # the brightness values should be def test_compute_posteriors(): @@ -113,7 +112,5 @@ def test_compute_posteriors(): if __name__ == "__main__": # test_brightness_calculation() - # test_read_input_with_brightness() - test_compute_posteriors() - - # Test commit + test_read_input_with_brightness() + # test_compute_posteriors() From 2ca8419ca66d652f88afdad0880b098e7517d622 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 20 Sep 2024 13:12:31 -0500 Subject: [PATCH 18/44] cleaned up test brightness --- tests/Visual.py | 50 ++++++++++++++++++++++++++++++++++++++++ tests/test_brightness.py | 11 +++++---- 2 files changed, 56 insertions(+), 5 deletions(-) create mode 100644 tests/Visual.py diff --git a/tests/Visual.py b/tests/Visual.py new file mode 100644 index 00000000..a56db702 --- /dev/null +++ b/tests/Visual.py @@ -0,0 +1,50 @@ +from orbitize import system, read_input, DATADIR +import os +import numpy as np +import matplotlib.pyplot as plt + + +num_secondary_bodies = 1 + +input_file = os.path.join(DATADIR, "GJ504.csv") +data_table = read_input.read_file(input_file) + +system_mass = 1.47 +plx = 24.30 + +test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) + +params_arr = np.array( + [ + 10.0, # sma + 0.9, # ecc + np.radians(30), # inc + np.radians(60), # aop + np.radians(120), # pan + 0.0, # tau + plx, + system_mass, + ] +) +epochs = np.linspace(0, 365 * 30, int(1e3)) +ra, dec, vz, brightness = test_system.compute_all_orbits(params_arr, epochs=epochs) + +fig, ax = plt.subplots(2, 1, figsize=(5, 10)) + +ax[0].scatter( + epochs, + brightness, + color=plt.cm.RdYlBu((epochs - epochs[0]) / (epochs[-1] - epochs[0])), +) + +ax[1].scatter( + ra[:, 1, :], + dec[:, 1, :], + color=plt.cm.RdYlBu((epochs - epochs[0]) / (epochs[-1] - epochs[0])), +) + +ax[1].scatter([0], [0], color="red") + +ax[1].axis("equal") + +plt.savefig("visual4farrah.png") \ No newline at end of file diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 9c5ca348..e7b094ee 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -29,23 +29,25 @@ def test_brightness_calculation(): params = np.array( [ 10.0, - 0.1, + 0.3, np.radians(89), np.radians(21), np.radians(31), - 0.0, # note: I didn't convert tau here, just picked random number + 0.0, 51.5, 1.75, ] ) ra, dec, vz, brightness = test_system.compute_all_orbits(params) + + # TODO (farrah): make plot of brightness vs time plt.figure() plt.scatter(times, brightness) - plt.xlabel("Time [dy]") - plt.ylabel("Brightness") + plt.xlabel("Time [dy]", fontsize=18) + plt.ylabel("Brightness", fontsize=18) plt.savefig("Test_brightness.png") @@ -62,7 +64,6 @@ def test_read_input_with_brightness(): times = data_table["epoch"].value brightness_values = data_table["brightness"].value - # Do we need the rest of this? since the values for time and brightness are given print("hello! :D ") From 30248b652289d385b15bdda1dfacfeb16b7ba8f6 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Fri, 20 Sep 2024 11:41:29 -0700 Subject: [PATCH 19/44] testing add change --- tests/test_brightness.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 464dec53..05a04702 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -8,6 +8,8 @@ import numpy as np import matplotlib.pyplot as plt +print('hello') + def test_brightness_calculation(): From aa16ef82f6a563ed6e045f28a6f392c39e5af891 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 27 Sep 2024 13:39:00 -0500 Subject: [PATCH 20/44] Added assert brightness test --- tests/test_brightness.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 58fba946..7d6af83d 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -64,6 +64,22 @@ def test_read_input_with_brightness(): # TODO (Farrah): add a test that asserts the brightness column of the data table is # what you expect (hint: check in the reflected_light_example.csv to see what # the brightness values should be +def test_assert_brightness(): + + num_secondary_bodies = 1 + + input_file = os.path.join(DATADIR, "reflected_light_example.csv") + + data_table = read_input.read_file(input_file) + + brightness_values = data_table["brightness"].value + assert 'brightness' in data_table.columns, "Brightness column does not exist" + print("The assert works!") + print (brightness_values) + x = "welcome" + + #if condition returns False, AssertionError is raised: + assert x != "hello", "This is meant to say welcome" def test_compute_posteriors(): @@ -112,5 +128,6 @@ def test_compute_posteriors(): if __name__ == "__main__": # test_brightness_calculation() - test_read_input_with_brightness() + #test_read_input_with_brightness() + test_assert_brightness() # test_compute_posteriors() From cd0c00a5386e9c7ca5ab0fc9ca2b8aa922f8480d Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 8 Oct 2024 16:13:07 -0700 Subject: [PATCH 21/44] mcmc running --- orbitize/system.py | 7 +++++++ tests/test_brightness.py | 11 +++++------ 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index c016b4ba..d2ed093c 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -98,10 +98,15 @@ def __init__( # List of index arrays corresponding to each rv for each body self.rv = [] + # index arrays corresponding to brightness for each body + self.brightness = [] + self.fit_astrometry = True radec_indices = np.where(self.data_table["quant_type"] == "radec") seppa_indices = np.where(self.data_table["quant_type"] == "seppa") + brightness_indices = np.where(~np.isnan(self.data_table["brightness"])) + if len(radec_indices[0]) == 0 and len(seppa_indices[0]) == 0: self.fit_astrometry = False rv_indices = np.where(self.data_table["quant_type"] == "rv") @@ -140,6 +145,7 @@ def __init__( np.intersect1d(self.body_indices[body_num], seppa_indices) ) self.rv.append(np.intersect1d(self.body_indices[body_num], rv_indices)) + self.brightness.append(np.intersect1d(self.body_indices[body_num], brightness_indices)) # we should track the influence of the planet(s) on each other/the star if: # we are not fitting massless planets and @@ -302,6 +308,7 @@ def __init__( self.param_idx = self.basis.param_idx + def save(self, hf): """ Saves the current object to an hdf5 file diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 36afe73e..7f938bdf 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -120,16 +120,15 @@ def test_compute_posteriors(): ax[1].axis('equal') plt.savefig('visual4farrah.png') - # model = test_system.compute_model(params_arr) + model = test_system.compute_model(params_arr) # print(model) - # test_mcmc = sampler.MCMC(test_system, 1, 50, num_threads=1) - - # test_mcmc.run_sampler(10) + test_mcmc = sampler.MCMC(test_system, 1, 50, num_threads=1) + test_mcmc.run_sampler(10) if __name__ == "__main__": # test_brightness_calculation() #test_read_input_with_brightness() - test_assert_brightness() - # test_compute_posteriors() + # test_assert_brightness() + test_compute_posteriors() From 1a069bb6facfa598668a468556bd8298573edae6 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Wed, 9 Oct 2024 13:01:55 -0500 Subject: [PATCH 22/44] Assert nan values --- tests/test_brightness.py | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 7d6af83d..a1b4ca12 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -6,6 +6,7 @@ from orbitize import DATADIR, read_input import os import numpy as np +import pandas as pd import matplotlib.pyplot as plt @@ -32,7 +33,7 @@ def test_brightness_calculation(): 0.3, np.radians(89), np.radians(21), - np.radians(31), + np.radians(70), 0.0, 51.5, 1.75, @@ -45,7 +46,7 @@ def test_brightness_calculation(): plt.scatter(times, brightness) plt.xlabel("Time [dy]", fontsize=18) plt.ylabel("Brightness", fontsize=18) - plt.savefig("Test_brightness.png") + plt.savefig("Test_brightness2.png") def test_read_input_with_brightness(): @@ -74,13 +75,20 @@ def test_assert_brightness(): brightness_values = data_table["brightness"].value assert 'brightness' in data_table.columns, "Brightness column does not exist" - print("The assert works!") - print (brightness_values) - x = "welcome" + #print("The assert works!") + #print (brightness_values) + - #if condition returns False, AssertionError is raised: - assert x != "hello", "This is meant to say welcome" +def test_assert_nan(): + num_secondary_bodies = 1 + + input_file = os.path.join(DATADIR, "reflected_light_example.csv") + data_table = read_input.read_file(input_file) + + brightness_values = data_table["brightness"].value + assert pd.isna(data_table["brightness"][4]), "This is supposed to say NAN" + print(brightness_values) def test_compute_posteriors(): @@ -127,7 +135,8 @@ def test_compute_posteriors(): if __name__ == "__main__": - # test_brightness_calculation() + #test_brightness_calculation() #test_read_input_with_brightness() - test_assert_brightness() + #test_assert_brightness() + test_assert_nan() # test_compute_posteriors() From 73000151f85ab3e1e4cbf2e87220e67d438367dc Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Thu, 10 Oct 2024 11:15:41 -0500 Subject: [PATCH 23/44] changed these to test 'test_brightness.py' --- orbitize/example_data/reflected_light_example.csv | 2 +- tests/Visual.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/orbitize/example_data/reflected_light_example.csv b/orbitize/example_data/reflected_light_example.csv index 5a17cc4f..8f0fcec9 100644 --- a/orbitize/example_data/reflected_light_example.csv +++ b/orbitize/example_data/reflected_light_example.csv @@ -3,7 +3,7 @@ epoch,object,raoff,decoff,raoff_err,decoff_err,brightness 57606,1,236.63,127.94,9.77,9.18,0.5 57645,1,234.52,123.39,1.79,1.03,0.5 57946,1,210.76,152.09,1.94,1.88,0.5 -58276,1,167.49,180.87,1.61,16.97,0.5 +58276,1,167.49,180.87,1.61,16.97, 58287,1,177.67,174.6,1.67,1.67,0.5 58365,1,165.7,185.33,3.28,3.66,0.5 58368,1,170.38,185.94,2.52,2.74,0.5 diff --git a/tests/Visual.py b/tests/Visual.py index a56db702..6e418eb9 100644 --- a/tests/Visual.py +++ b/tests/Visual.py @@ -18,9 +18,9 @@ [ 10.0, # sma 0.9, # ecc - np.radians(30), # inc + np.radians(90), # inc np.radians(60), # aop - np.radians(120), # pan + np.radians(90), # pan 0.0, # tau plx, system_mass, From 683094aa544a42b46184ca0e285e8fe792c25230 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 12:20:41 -0500 Subject: [PATCH 24/44] Test orbital data to run mcmc --- test_orbital_data.csv | 51 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 test_orbital_data.csv diff --git a/test_orbital_data.csv b/test_orbital_data.csv new file mode 100644 index 00000000..71bc14b8 --- /dev/null +++ b/test_orbital_data.csv @@ -0,0 +1,51 @@ +raoff,decoff,brightness,epoch +0.811231477773281,0.7594337387611091,0.09358490399640425,0.0 +0.8086976534668815,0.7884882174763149,0.09862733963920862,51.02040816326531 +0.8018051182543783,0.8135752998448389,0.10499582237388896,102.04081632653062 +0.7904122778426745,0.8345316311433326,0.11302378885599486,153.06122448979593 +0.7743165421613408,0.8511262472797302,0.12318805546636806,204.08163265306123 +0.753249548414535,0.8630503873758559,0.1361800248346027,255.10204081632654 +0.7268712100315663,0.8699046823746542,0.1530226526335016,306.12244897959187 +0.694762808612201,0.871183254587578,0.17526767302647034,357.14285714285717 +0.6564196199390899,0.8662543021471459,0.2053375622554458,408.16326530612247 +0.6112440678303533,0.8543369207025031,0.2471342724034155,459.18367346938777 +0.5585412933875817,0.8344743818340995,0.30714302736709537,510.2040816326531 +0.49752057854752274,0.8055051180185213,0.39642605066463393,561.2244897959184 +0.42730869340142713,0.7660347500990113,0.5339941531274536,612.2448979591837 +0.3469855704957162,0.7144164850051454,0.7510921949147653,663.265306122449 +0.25565954973242283,0.6487544841765712,1.089659065489047,714.2857142857143 +0.15260947588849744,0.5669573266103906,1.5668452007999316,765.3061224489796 +0.037533685832437026,0.4668885487044668,2.062626571880702,816.3265306122449 +-0.08904312701616794,0.34668838046440426,2.2566005333932724,867.3469387755102 +-0.22515974344986556,0.2053664293775624,1.9742762557365308,918.3673469387755 +-0.36661560366594714,0.04375763921362707,1.4724201265196528,969.3877551020408 +-0.5061474469043882,-0.13417842573759436,1.0378267478370446,1020.4081632653061 +-0.6330566925723007,-0.32005017911521655,0.7457820123781923,1071.4285714285716 +-0.7340642827698997,-0.5005631649102175,0.5692403465256193,1122.4489795918369 +-0.7960233988208258,-0.6590354524244593,0.3677023181237994,1173.4693877551022 +-0.8100093753970008,-0.779341500288983,0.310805460038203,1224.4897959183675 +-0.7747106962319243,-0.8508197707379554,0.41184810403842464,1275.5102040816328 +-0.6968201319288353,-0.8712731010397949,0.4267037403343691,1326.530612244898 +-0.5881886796558359,-0.8463291454704698,0.45692118233449847,1377.5510204081634 +-0.46179505907320606,-0.7861524685145167,0.4835731029852539,1428.5714285714287 +-0.32874319217715686,-0.7018568896193422,0.47696040679432183,1479.591836734694 +-0.19704231837725736,-0.603169261080095,0.4173135946760028,1530.6122448979593 +-0.0717085052125464,-0.49753210156043054,0.32665907139175554,1581.6326530612246 +0.04455725697773338,-0.39013739669169883,0.24572306795533427,1632.6530612244899 +0.15061964729048624,-0.2843546293265123,0.189864891641152,1683.6734693877552 +0.246302320136533,-0.18222582266793214,0.1541293228877271,1734.6938775510205 +0.3319636092307561,-0.08488826595156958,0.1309447698576423,1785.7142857142858 +0.4082247422781462,0.007108706138521259,0.11532666162626576,1836.734693877551 +0.47580638683140664,0.09358245909783232,0.10445257538616691,1887.7551020408164 +0.5354357374912007,0.17456889152858981,0.09672357063471888,1938.7755102040817 +0.5877970769939762,0.25022587899384413,0.09119647427470219,1989.795918367347 +0.6335080358663165,0.3207715550593558,0.08729085845318117,2040.8163265306123 +0.673110446052766,0.3864454566166559,0.08463812457943871,2091.8367346938776 +0.707069071798058,0.44748437788656903,0.0830008016213394,2142.857142857143 +0.7357742511645264,0.5041075168486121,0.08222762558027762,2193.877551020408 +0.75954615418438,0.5565073629841016,0.08222788911616931,2244.8979591836737 +0.7786393591852832,0.6048440073048987,0.08295689503656585,2295.918367346939 +0.7932470300017178,0.6492413535736495,0.0844083805697474,2346.9387755102043 +0.8035043081850315,0.6897842167972191,0.08661187268900492,2397.9591836734694 +0.8094907181927744,0.7265156142343624,0.08963415615760459,2448.979591836735 +0.811231477773281,0.7594337387611084,0.09358490399640412,2500.0 From 1a6f24e8f85e577a43741fa3be2766b26222cb95 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 12:27:04 -0500 Subject: [PATCH 25/44] Moved test orbital data to example data file --- orbitize/example_data/test_orbital_data.csv | 51 +++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 orbitize/example_data/test_orbital_data.csv diff --git a/orbitize/example_data/test_orbital_data.csv b/orbitize/example_data/test_orbital_data.csv new file mode 100644 index 00000000..71bc14b8 --- /dev/null +++ b/orbitize/example_data/test_orbital_data.csv @@ -0,0 +1,51 @@ +raoff,decoff,brightness,epoch +0.811231477773281,0.7594337387611091,0.09358490399640425,0.0 +0.8086976534668815,0.7884882174763149,0.09862733963920862,51.02040816326531 +0.8018051182543783,0.8135752998448389,0.10499582237388896,102.04081632653062 +0.7904122778426745,0.8345316311433326,0.11302378885599486,153.06122448979593 +0.7743165421613408,0.8511262472797302,0.12318805546636806,204.08163265306123 +0.753249548414535,0.8630503873758559,0.1361800248346027,255.10204081632654 +0.7268712100315663,0.8699046823746542,0.1530226526335016,306.12244897959187 +0.694762808612201,0.871183254587578,0.17526767302647034,357.14285714285717 +0.6564196199390899,0.8662543021471459,0.2053375622554458,408.16326530612247 +0.6112440678303533,0.8543369207025031,0.2471342724034155,459.18367346938777 +0.5585412933875817,0.8344743818340995,0.30714302736709537,510.2040816326531 +0.49752057854752274,0.8055051180185213,0.39642605066463393,561.2244897959184 +0.42730869340142713,0.7660347500990113,0.5339941531274536,612.2448979591837 +0.3469855704957162,0.7144164850051454,0.7510921949147653,663.265306122449 +0.25565954973242283,0.6487544841765712,1.089659065489047,714.2857142857143 +0.15260947588849744,0.5669573266103906,1.5668452007999316,765.3061224489796 +0.037533685832437026,0.4668885487044668,2.062626571880702,816.3265306122449 +-0.08904312701616794,0.34668838046440426,2.2566005333932724,867.3469387755102 +-0.22515974344986556,0.2053664293775624,1.9742762557365308,918.3673469387755 +-0.36661560366594714,0.04375763921362707,1.4724201265196528,969.3877551020408 +-0.5061474469043882,-0.13417842573759436,1.0378267478370446,1020.4081632653061 +-0.6330566925723007,-0.32005017911521655,0.7457820123781923,1071.4285714285716 +-0.7340642827698997,-0.5005631649102175,0.5692403465256193,1122.4489795918369 +-0.7960233988208258,-0.6590354524244593,0.3677023181237994,1173.4693877551022 +-0.8100093753970008,-0.779341500288983,0.310805460038203,1224.4897959183675 +-0.7747106962319243,-0.8508197707379554,0.41184810403842464,1275.5102040816328 +-0.6968201319288353,-0.8712731010397949,0.4267037403343691,1326.530612244898 +-0.5881886796558359,-0.8463291454704698,0.45692118233449847,1377.5510204081634 +-0.46179505907320606,-0.7861524685145167,0.4835731029852539,1428.5714285714287 +-0.32874319217715686,-0.7018568896193422,0.47696040679432183,1479.591836734694 +-0.19704231837725736,-0.603169261080095,0.4173135946760028,1530.6122448979593 +-0.0717085052125464,-0.49753210156043054,0.32665907139175554,1581.6326530612246 +0.04455725697773338,-0.39013739669169883,0.24572306795533427,1632.6530612244899 +0.15061964729048624,-0.2843546293265123,0.189864891641152,1683.6734693877552 +0.246302320136533,-0.18222582266793214,0.1541293228877271,1734.6938775510205 +0.3319636092307561,-0.08488826595156958,0.1309447698576423,1785.7142857142858 +0.4082247422781462,0.007108706138521259,0.11532666162626576,1836.734693877551 +0.47580638683140664,0.09358245909783232,0.10445257538616691,1887.7551020408164 +0.5354357374912007,0.17456889152858981,0.09672357063471888,1938.7755102040817 +0.5877970769939762,0.25022587899384413,0.09119647427470219,1989.795918367347 +0.6335080358663165,0.3207715550593558,0.08729085845318117,2040.8163265306123 +0.673110446052766,0.3864454566166559,0.08463812457943871,2091.8367346938776 +0.707069071798058,0.44748437788656903,0.0830008016213394,2142.857142857143 +0.7357742511645264,0.5041075168486121,0.08222762558027762,2193.877551020408 +0.75954615418438,0.5565073629841016,0.08222788911616931,2244.8979591836737 +0.7786393591852832,0.6048440073048987,0.08295689503656585,2295.918367346939 +0.7932470300017178,0.6492413535736495,0.0844083805697474,2346.9387755102043 +0.8035043081850315,0.6897842167972191,0.08661187268900492,2397.9591836734694 +0.8094907181927744,0.7265156142343624,0.08963415615760459,2448.979591836735 +0.811231477773281,0.7594337387611084,0.09358490399640412,2500.0 From 32341952e7af5432ee2003d98a1cc6487d593744 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 12:47:36 -0500 Subject: [PATCH 26/44] Added ID to test table --- .../example_data/orbital_data_with_id.csv | 51 +++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 orbitize/example_data/orbital_data_with_id.csv diff --git a/orbitize/example_data/orbital_data_with_id.csv b/orbitize/example_data/orbital_data_with_id.csv new file mode 100644 index 00000000..e3b5024c --- /dev/null +++ b/orbitize/example_data/orbital_data_with_id.csv @@ -0,0 +1,51 @@ +object_id,ra,dec,brightness,epoch +1,0.811231477773281,0.7594337387611091,0.09358490399640425,0.0 +2,0.8086976534668815,0.7884882174763149,0.09862733963920862,51.02040816326531 +3,0.8018051182543783,0.8135752998448389,0.10499582237388896,102.04081632653062 +4,0.7904122778426745,0.8345316311433326,0.11302378885599486,153.06122448979593 +5,0.7743165421613408,0.8511262472797302,0.12318805546636806,204.08163265306123 +6,0.753249548414535,0.8630503873758559,0.1361800248346027,255.10204081632654 +7,0.7268712100315663,0.8699046823746542,0.1530226526335016,306.12244897959187 +8,0.694762808612201,0.871183254587578,0.17526767302647034,357.14285714285717 +9,0.6564196199390899,0.8662543021471459,0.2053375622554458,408.16326530612247 +10,0.6112440678303533,0.8543369207025031,0.2471342724034155,459.18367346938777 +11,0.5585412933875817,0.8344743818340995,0.30714302736709537,510.2040816326531 +12,0.49752057854752274,0.8055051180185213,0.39642605066463393,561.2244897959184 +13,0.42730869340142713,0.7660347500990113,0.5339941531274536,612.2448979591837 +14,0.3469855704957162,0.7144164850051454,0.7510921949147653,663.265306122449 +15,0.25565954973242283,0.6487544841765712,1.089659065489047,714.2857142857143 +16,0.15260947588849744,0.5669573266103906,1.5668452007999316,765.3061224489796 +17,0.037533685832437026,0.4668885487044668,2.062626571880702,816.3265306122449 +18,-0.08904312701616794,0.34668838046440426,2.2566005333932724,867.3469387755102 +19,-0.22515974344986556,0.2053664293775624,1.9742762557365308,918.3673469387755 +20,-0.36661560366594714,0.04375763921362707,1.4724201265196528,969.3877551020408 +21,-0.5061474469043882,-0.13417842573759436,1.0378267478370446,1020.4081632653061 +22,-0.6330566925723007,-0.32005017911521655,0.7457820123781923,1071.4285714285716 +23,-0.7340642827698997,-0.5005631649102175,0.5692403465256193,1122.4489795918369 +24,-0.7960233988208258,-0.6590354524244593,0.3677023181237994,1173.4693877551022 +25,-0.8100093753970008,-0.779341500288983,0.310805460038203,1224.4897959183675 +26,-0.7747106962319243,-0.8508197707379554,0.41184810403842464,1275.5102040816328 +27,-0.6968201319288353,-0.8712731010397949,0.4267037403343691,1326.530612244898 +28,-0.5881886796558359,-0.8463291454704698,0.45692118233449847,1377.5510204081634 +29,-0.46179505907320606,-0.7861524685145167,0.4835731029852539,1428.5714285714287 +30,-0.32874319217715686,-0.7018568896193422,0.47696040679432183,1479.591836734694 +31,-0.19704231837725736,-0.603169261080095,0.4173135946760028,1530.6122448979593 +32,-0.0717085052125464,-0.49753210156043054,0.32665907139175554,1581.6326530612246 +33,0.04455725697773338,-0.39013739669169883,0.24572306795533427,1632.6530612244899 +34,0.15061964729048624,-0.2843546293265123,0.189864891641152,1683.6734693877552 +35,0.246302320136533,-0.18222582266793214,0.1541293228877271,1734.6938775510205 +36,0.3319636092307561,-0.08488826595156958,0.1309447698576423,1785.7142857142858 +37,0.4082247422781462,0.007108706138521259,0.11532666162626576,1836.734693877551 +38,0.47580638683140664,0.09358245909783232,0.10445257538616691,1887.7551020408164 +39,0.5354357374912007,0.17456889152858981,0.09672357063471888,1938.7755102040817 +40,0.5877970769939762,0.25022587899384413,0.09119647427470219,1989.795918367347 +41,0.6335080358663165,0.3207715550593558,0.08729085845318117,2040.8163265306123 +42,0.673110446052766,0.3864454566166559,0.08463812457943871,2091.8367346938776 +43,0.707069071798058,0.44748437788656903,0.0830008016213394,2142.857142857143 +44,0.7357742511645264,0.5041075168486121,0.08222762558027762,2193.877551020408 +45,0.75954615418438,0.5565073629841016,0.08222788911616931,2244.8979591836737 +46,0.7786393591852832,0.6048440073048987,0.08295689503656585,2295.918367346939 +47,0.7932470300017178,0.6492413535736495,0.0844083805697474,2346.9387755102043 +48,0.8035043081850315,0.6897842167972191,0.08661187268900492,2397.9591836734694 +49,0.8094907181927744,0.7265156142343624,0.08963415615760459,2448.979591836735 +50,0.811231477773281,0.7594337387611084,0.09358490399640412,2500.0 From ac4ac17d793a5b39ac6ee53cbf6c1c48312e2dac Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 13:07:10 -0500 Subject: [PATCH 27/44] another change that has columns in right order...lol --- .../example_data/orbital_data_with_id.csv | 102 +++++++++--------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/orbitize/example_data/orbital_data_with_id.csv b/orbitize/example_data/orbital_data_with_id.csv index e3b5024c..1c8edd70 100644 --- a/orbitize/example_data/orbital_data_with_id.csv +++ b/orbitize/example_data/orbital_data_with_id.csv @@ -1,51 +1,51 @@ -object_id,ra,dec,brightness,epoch -1,0.811231477773281,0.7594337387611091,0.09358490399640425,0.0 -2,0.8086976534668815,0.7884882174763149,0.09862733963920862,51.02040816326531 -3,0.8018051182543783,0.8135752998448389,0.10499582237388896,102.04081632653062 -4,0.7904122778426745,0.8345316311433326,0.11302378885599486,153.06122448979593 -5,0.7743165421613408,0.8511262472797302,0.12318805546636806,204.08163265306123 -6,0.753249548414535,0.8630503873758559,0.1361800248346027,255.10204081632654 -7,0.7268712100315663,0.8699046823746542,0.1530226526335016,306.12244897959187 -8,0.694762808612201,0.871183254587578,0.17526767302647034,357.14285714285717 -9,0.6564196199390899,0.8662543021471459,0.2053375622554458,408.16326530612247 -10,0.6112440678303533,0.8543369207025031,0.2471342724034155,459.18367346938777 -11,0.5585412933875817,0.8344743818340995,0.30714302736709537,510.2040816326531 -12,0.49752057854752274,0.8055051180185213,0.39642605066463393,561.2244897959184 -13,0.42730869340142713,0.7660347500990113,0.5339941531274536,612.2448979591837 -14,0.3469855704957162,0.7144164850051454,0.7510921949147653,663.265306122449 -15,0.25565954973242283,0.6487544841765712,1.089659065489047,714.2857142857143 -16,0.15260947588849744,0.5669573266103906,1.5668452007999316,765.3061224489796 -17,0.037533685832437026,0.4668885487044668,2.062626571880702,816.3265306122449 -18,-0.08904312701616794,0.34668838046440426,2.2566005333932724,867.3469387755102 -19,-0.22515974344986556,0.2053664293775624,1.9742762557365308,918.3673469387755 -20,-0.36661560366594714,0.04375763921362707,1.4724201265196528,969.3877551020408 -21,-0.5061474469043882,-0.13417842573759436,1.0378267478370446,1020.4081632653061 -22,-0.6330566925723007,-0.32005017911521655,0.7457820123781923,1071.4285714285716 -23,-0.7340642827698997,-0.5005631649102175,0.5692403465256193,1122.4489795918369 -24,-0.7960233988208258,-0.6590354524244593,0.3677023181237994,1173.4693877551022 -25,-0.8100093753970008,-0.779341500288983,0.310805460038203,1224.4897959183675 -26,-0.7747106962319243,-0.8508197707379554,0.41184810403842464,1275.5102040816328 -27,-0.6968201319288353,-0.8712731010397949,0.4267037403343691,1326.530612244898 -28,-0.5881886796558359,-0.8463291454704698,0.45692118233449847,1377.5510204081634 -29,-0.46179505907320606,-0.7861524685145167,0.4835731029852539,1428.5714285714287 -30,-0.32874319217715686,-0.7018568896193422,0.47696040679432183,1479.591836734694 -31,-0.19704231837725736,-0.603169261080095,0.4173135946760028,1530.6122448979593 -32,-0.0717085052125464,-0.49753210156043054,0.32665907139175554,1581.6326530612246 -33,0.04455725697773338,-0.39013739669169883,0.24572306795533427,1632.6530612244899 -34,0.15061964729048624,-0.2843546293265123,0.189864891641152,1683.6734693877552 -35,0.246302320136533,-0.18222582266793214,0.1541293228877271,1734.6938775510205 -36,0.3319636092307561,-0.08488826595156958,0.1309447698576423,1785.7142857142858 -37,0.4082247422781462,0.007108706138521259,0.11532666162626576,1836.734693877551 -38,0.47580638683140664,0.09358245909783232,0.10445257538616691,1887.7551020408164 -39,0.5354357374912007,0.17456889152858981,0.09672357063471888,1938.7755102040817 -40,0.5877970769939762,0.25022587899384413,0.09119647427470219,1989.795918367347 -41,0.6335080358663165,0.3207715550593558,0.08729085845318117,2040.8163265306123 -42,0.673110446052766,0.3864454566166559,0.08463812457943871,2091.8367346938776 -43,0.707069071798058,0.44748437788656903,0.0830008016213394,2142.857142857143 -44,0.7357742511645264,0.5041075168486121,0.08222762558027762,2193.877551020408 -45,0.75954615418438,0.5565073629841016,0.08222788911616931,2244.8979591836737 -46,0.7786393591852832,0.6048440073048987,0.08295689503656585,2295.918367346939 -47,0.7932470300017178,0.6492413535736495,0.0844083805697474,2346.9387755102043 -48,0.8035043081850315,0.6897842167972191,0.08661187268900492,2397.9591836734694 -49,0.8094907181927744,0.7265156142343624,0.08963415615760459,2448.979591836735 -50,0.811231477773281,0.7594337387611084,0.09358490399640412,2500.0 +epoch,object_id,raoff,decoff,brightness +0.0,1,0.811231477773281,0.7594337387611091,0.09358490399640425 +51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.09862733963920862 +102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.10499582237388896 +153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.11302378885599486 +204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.12318805546636806 +255.10204081632654,1,0.753249548414535,0.8630503873758559,0.1361800248346027 +306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.1530226526335016 +357.14285714285717,1,0.694762808612201,0.871183254587578,0.17526767302647034 +408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.2053375622554458 +459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.2471342724034155 +510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.30714302736709537 +561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.39642605066463393 +612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.5339941531274536 +663.265306122449,1,0.3469855704957162,0.7144164850051454,0.7510921949147653 +714.2857142857143,1,0.25565954973242283,0.6487544841765712,1.089659065489047 +765.3061224489796,1,0.15260947588849744,0.5669573266103906,1.5668452007999316 +816.3265306122449,1,0.037533685832437026,0.4668885487044668,2.062626571880702 +867.3469387755102,1,-0.08904312701616794,0.34668838046440426,2.2566005333932724 +918.3673469387755,1,-0.22515974344986556,0.2053664293775624,1.9742762557365308 +969.3877551020408,1,-0.36661560366594714,0.04375763921362707,1.4724201265196528 +1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,1.0378267478370446 +1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.7457820123781923 +1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.5692403465256193 +1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.3677023181237994 +1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.310805460038203 +1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.41184810403842464 +1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.4267037403343691 +1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.45692118233449847 +1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.4835731029852539 +1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.47696040679432183 +1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.4173135946760028 +1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.32665907139175554 +1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.24572306795533427 +1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.189864891641152 +1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.1541293228877271 +1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.1309447698576423 +1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.11532666162626576 +1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.10445257538616691 +1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.09672357063471888 +1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.09119647427470219 +2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.08729085845318117 +2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.08463812457943871 +2142.857142857143,1,0.707069071798058,0.44748437788656903,0.0830008016213394 +2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.08222762558027762 +2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.08222788911616931 +2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.08295689503656585 +2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.0844083805697474 +2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.08661187268900492 +2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.08963415615760459 +2500.0,1,0.811231477773281,0.7594337387611084,0.09358490399640412 From 3147d1b732d24d201c2d8ca4041b4c821addf03e Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 13:09:19 -0500 Subject: [PATCH 28/44] object_id changed to object to fit exception! --- orbitize/example_data/orbital_data_with_id.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orbitize/example_data/orbital_data_with_id.csv b/orbitize/example_data/orbital_data_with_id.csv index 1c8edd70..47efabda 100644 --- a/orbitize/example_data/orbital_data_with_id.csv +++ b/orbitize/example_data/orbital_data_with_id.csv @@ -1,4 +1,4 @@ -epoch,object_id,raoff,decoff,brightness +epoch,object,raoff,decoff,brightness 0.0,1,0.811231477773281,0.7594337387611091,0.09358490399640425 51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.09862733963920862 102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.10499582237388896 From f40b8db4b4527243667f25e52e9b7cd4c8554371 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 13:13:38 -0500 Subject: [PATCH 29/44] added + read in a simulated data file to compute posteriors test --- orbitize/example_data/test_orbital_data.csv | 51 --------------------- test_orbital_data.csv | 51 --------------------- tests/test_brightness.py | 4 +- 3 files changed, 2 insertions(+), 104 deletions(-) delete mode 100644 orbitize/example_data/test_orbital_data.csv delete mode 100644 test_orbital_data.csv diff --git a/orbitize/example_data/test_orbital_data.csv b/orbitize/example_data/test_orbital_data.csv deleted file mode 100644 index 71bc14b8..00000000 --- a/orbitize/example_data/test_orbital_data.csv +++ /dev/null @@ -1,51 +0,0 @@ -raoff,decoff,brightness,epoch -0.811231477773281,0.7594337387611091,0.09358490399640425,0.0 -0.8086976534668815,0.7884882174763149,0.09862733963920862,51.02040816326531 -0.8018051182543783,0.8135752998448389,0.10499582237388896,102.04081632653062 -0.7904122778426745,0.8345316311433326,0.11302378885599486,153.06122448979593 -0.7743165421613408,0.8511262472797302,0.12318805546636806,204.08163265306123 -0.753249548414535,0.8630503873758559,0.1361800248346027,255.10204081632654 -0.7268712100315663,0.8699046823746542,0.1530226526335016,306.12244897959187 -0.694762808612201,0.871183254587578,0.17526767302647034,357.14285714285717 -0.6564196199390899,0.8662543021471459,0.2053375622554458,408.16326530612247 -0.6112440678303533,0.8543369207025031,0.2471342724034155,459.18367346938777 -0.5585412933875817,0.8344743818340995,0.30714302736709537,510.2040816326531 -0.49752057854752274,0.8055051180185213,0.39642605066463393,561.2244897959184 -0.42730869340142713,0.7660347500990113,0.5339941531274536,612.2448979591837 -0.3469855704957162,0.7144164850051454,0.7510921949147653,663.265306122449 -0.25565954973242283,0.6487544841765712,1.089659065489047,714.2857142857143 -0.15260947588849744,0.5669573266103906,1.5668452007999316,765.3061224489796 -0.037533685832437026,0.4668885487044668,2.062626571880702,816.3265306122449 --0.08904312701616794,0.34668838046440426,2.2566005333932724,867.3469387755102 --0.22515974344986556,0.2053664293775624,1.9742762557365308,918.3673469387755 --0.36661560366594714,0.04375763921362707,1.4724201265196528,969.3877551020408 --0.5061474469043882,-0.13417842573759436,1.0378267478370446,1020.4081632653061 --0.6330566925723007,-0.32005017911521655,0.7457820123781923,1071.4285714285716 --0.7340642827698997,-0.5005631649102175,0.5692403465256193,1122.4489795918369 --0.7960233988208258,-0.6590354524244593,0.3677023181237994,1173.4693877551022 --0.8100093753970008,-0.779341500288983,0.310805460038203,1224.4897959183675 --0.7747106962319243,-0.8508197707379554,0.41184810403842464,1275.5102040816328 --0.6968201319288353,-0.8712731010397949,0.4267037403343691,1326.530612244898 --0.5881886796558359,-0.8463291454704698,0.45692118233449847,1377.5510204081634 --0.46179505907320606,-0.7861524685145167,0.4835731029852539,1428.5714285714287 --0.32874319217715686,-0.7018568896193422,0.47696040679432183,1479.591836734694 --0.19704231837725736,-0.603169261080095,0.4173135946760028,1530.6122448979593 --0.0717085052125464,-0.49753210156043054,0.32665907139175554,1581.6326530612246 -0.04455725697773338,-0.39013739669169883,0.24572306795533427,1632.6530612244899 -0.15061964729048624,-0.2843546293265123,0.189864891641152,1683.6734693877552 -0.246302320136533,-0.18222582266793214,0.1541293228877271,1734.6938775510205 -0.3319636092307561,-0.08488826595156958,0.1309447698576423,1785.7142857142858 -0.4082247422781462,0.007108706138521259,0.11532666162626576,1836.734693877551 -0.47580638683140664,0.09358245909783232,0.10445257538616691,1887.7551020408164 -0.5354357374912007,0.17456889152858981,0.09672357063471888,1938.7755102040817 -0.5877970769939762,0.25022587899384413,0.09119647427470219,1989.795918367347 -0.6335080358663165,0.3207715550593558,0.08729085845318117,2040.8163265306123 -0.673110446052766,0.3864454566166559,0.08463812457943871,2091.8367346938776 -0.707069071798058,0.44748437788656903,0.0830008016213394,2142.857142857143 -0.7357742511645264,0.5041075168486121,0.08222762558027762,2193.877551020408 -0.75954615418438,0.5565073629841016,0.08222788911616931,2244.8979591836737 -0.7786393591852832,0.6048440073048987,0.08295689503656585,2295.918367346939 -0.7932470300017178,0.6492413535736495,0.0844083805697474,2346.9387755102043 -0.8035043081850315,0.6897842167972191,0.08661187268900492,2397.9591836734694 -0.8094907181927744,0.7265156142343624,0.08963415615760459,2448.979591836735 -0.811231477773281,0.7594337387611084,0.09358490399640412,2500.0 diff --git a/test_orbital_data.csv b/test_orbital_data.csv deleted file mode 100644 index 71bc14b8..00000000 --- a/test_orbital_data.csv +++ /dev/null @@ -1,51 +0,0 @@ -raoff,decoff,brightness,epoch -0.811231477773281,0.7594337387611091,0.09358490399640425,0.0 -0.8086976534668815,0.7884882174763149,0.09862733963920862,51.02040816326531 -0.8018051182543783,0.8135752998448389,0.10499582237388896,102.04081632653062 -0.7904122778426745,0.8345316311433326,0.11302378885599486,153.06122448979593 -0.7743165421613408,0.8511262472797302,0.12318805546636806,204.08163265306123 -0.753249548414535,0.8630503873758559,0.1361800248346027,255.10204081632654 -0.7268712100315663,0.8699046823746542,0.1530226526335016,306.12244897959187 -0.694762808612201,0.871183254587578,0.17526767302647034,357.14285714285717 -0.6564196199390899,0.8662543021471459,0.2053375622554458,408.16326530612247 -0.6112440678303533,0.8543369207025031,0.2471342724034155,459.18367346938777 -0.5585412933875817,0.8344743818340995,0.30714302736709537,510.2040816326531 -0.49752057854752274,0.8055051180185213,0.39642605066463393,561.2244897959184 -0.42730869340142713,0.7660347500990113,0.5339941531274536,612.2448979591837 -0.3469855704957162,0.7144164850051454,0.7510921949147653,663.265306122449 -0.25565954973242283,0.6487544841765712,1.089659065489047,714.2857142857143 -0.15260947588849744,0.5669573266103906,1.5668452007999316,765.3061224489796 -0.037533685832437026,0.4668885487044668,2.062626571880702,816.3265306122449 --0.08904312701616794,0.34668838046440426,2.2566005333932724,867.3469387755102 --0.22515974344986556,0.2053664293775624,1.9742762557365308,918.3673469387755 --0.36661560366594714,0.04375763921362707,1.4724201265196528,969.3877551020408 --0.5061474469043882,-0.13417842573759436,1.0378267478370446,1020.4081632653061 --0.6330566925723007,-0.32005017911521655,0.7457820123781923,1071.4285714285716 --0.7340642827698997,-0.5005631649102175,0.5692403465256193,1122.4489795918369 --0.7960233988208258,-0.6590354524244593,0.3677023181237994,1173.4693877551022 --0.8100093753970008,-0.779341500288983,0.310805460038203,1224.4897959183675 --0.7747106962319243,-0.8508197707379554,0.41184810403842464,1275.5102040816328 --0.6968201319288353,-0.8712731010397949,0.4267037403343691,1326.530612244898 --0.5881886796558359,-0.8463291454704698,0.45692118233449847,1377.5510204081634 --0.46179505907320606,-0.7861524685145167,0.4835731029852539,1428.5714285714287 --0.32874319217715686,-0.7018568896193422,0.47696040679432183,1479.591836734694 --0.19704231837725736,-0.603169261080095,0.4173135946760028,1530.6122448979593 --0.0717085052125464,-0.49753210156043054,0.32665907139175554,1581.6326530612246 -0.04455725697773338,-0.39013739669169883,0.24572306795533427,1632.6530612244899 -0.15061964729048624,-0.2843546293265123,0.189864891641152,1683.6734693877552 -0.246302320136533,-0.18222582266793214,0.1541293228877271,1734.6938775510205 -0.3319636092307561,-0.08488826595156958,0.1309447698576423,1785.7142857142858 -0.4082247422781462,0.007108706138521259,0.11532666162626576,1836.734693877551 -0.47580638683140664,0.09358245909783232,0.10445257538616691,1887.7551020408164 -0.5354357374912007,0.17456889152858981,0.09672357063471888,1938.7755102040817 -0.5877970769939762,0.25022587899384413,0.09119647427470219,1989.795918367347 -0.6335080358663165,0.3207715550593558,0.08729085845318117,2040.8163265306123 -0.673110446052766,0.3864454566166559,0.08463812457943871,2091.8367346938776 -0.707069071798058,0.44748437788656903,0.0830008016213394,2142.857142857143 -0.7357742511645264,0.5041075168486121,0.08222762558027762,2193.877551020408 -0.75954615418438,0.5565073629841016,0.08222788911616931,2244.8979591836737 -0.7786393591852832,0.6048440073048987,0.08295689503656585,2295.918367346939 -0.7932470300017178,0.6492413535736495,0.0844083805697474,2346.9387755102043 -0.8035043081850315,0.6897842167972191,0.08661187268900492,2397.9591836734694 -0.8094907181927744,0.7265156142343624,0.08963415615760459,2448.979591836735 -0.811231477773281,0.7594337387611084,0.09358490399640412,2500.0 diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 5073a55a..60091cd1 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -55,7 +55,7 @@ def test_read_input_with_brightness(): num_secondary_bodies = 1 - input_file = os.path.join(DATADIR, "reflected_light_example.csv") + input_file = os.path.join(DATADIR, "test_orbital_data.csv") data_table = read_input.read_file(input_file) @@ -96,7 +96,7 @@ def test_compute_posteriors(): num_secondary_bodies = 1 - input_file = os.path.join(DATADIR, "GJ504.csv") + input_file = os.path.join(DATADIR, "orbital_data_with_id.csv") data_table = read_input.read_file(input_file) system_mass = 1.47 From 062b1f38578695cda0ab8c5863b8409f30071320 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 13:24:21 -0500 Subject: [PATCH 30/44] 100th attempt at fixing the csv file: adding ra and dec error --- .../example_data/orbital_data_with_id.csv | 102 +++++++++--------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/orbitize/example_data/orbital_data_with_id.csv b/orbitize/example_data/orbital_data_with_id.csv index 47efabda..5560a0f3 100644 --- a/orbitize/example_data/orbital_data_with_id.csv +++ b/orbitize/example_data/orbital_data_with_id.csv @@ -1,51 +1,51 @@ -epoch,object,raoff,decoff,brightness -0.0,1,0.811231477773281,0.7594337387611091,0.09358490399640425 -51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.09862733963920862 -102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.10499582237388896 -153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.11302378885599486 -204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.12318805546636806 -255.10204081632654,1,0.753249548414535,0.8630503873758559,0.1361800248346027 -306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.1530226526335016 -357.14285714285717,1,0.694762808612201,0.871183254587578,0.17526767302647034 -408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.2053375622554458 -459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.2471342724034155 -510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.30714302736709537 -561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.39642605066463393 -612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.5339941531274536 -663.265306122449,1,0.3469855704957162,0.7144164850051454,0.7510921949147653 -714.2857142857143,1,0.25565954973242283,0.6487544841765712,1.089659065489047 -765.3061224489796,1,0.15260947588849744,0.5669573266103906,1.5668452007999316 -816.3265306122449,1,0.037533685832437026,0.4668885487044668,2.062626571880702 -867.3469387755102,1,-0.08904312701616794,0.34668838046440426,2.2566005333932724 -918.3673469387755,1,-0.22515974344986556,0.2053664293775624,1.9742762557365308 -969.3877551020408,1,-0.36661560366594714,0.04375763921362707,1.4724201265196528 -1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,1.0378267478370446 -1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.7457820123781923 -1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.5692403465256193 -1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.3677023181237994 -1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.310805460038203 -1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.41184810403842464 -1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.4267037403343691 -1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.45692118233449847 -1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.4835731029852539 -1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.47696040679432183 -1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.4173135946760028 -1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.32665907139175554 -1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.24572306795533427 -1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.189864891641152 -1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.1541293228877271 -1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.1309447698576423 -1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.11532666162626576 -1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.10445257538616691 -1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.09672357063471888 -1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.09119647427470219 -2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.08729085845318117 -2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.08463812457943871 -2142.857142857143,1,0.707069071798058,0.44748437788656903,0.0830008016213394 -2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.08222762558027762 -2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.08222788911616931 -2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.08295689503656585 -2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.0844083805697474 -2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.08661187268900492 -2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.08963415615760459 -2500.0,1,0.811231477773281,0.7594337387611084,0.09358490399640412 +epoch,object,raoff,decoff,raoff_err,decoff_err,brightness +0.0,1,0.811231477773281,0.7594337387611091,nan,nan,0.09358490399640425 +51.02040816326531,1,0.8086976534668815,0.7884882174763149,nan,nan,0.09862733963920862 +102.04081632653062,1,0.8018051182543783,0.8135752998448389,nan,nan,0.10499582237388896 +153.06122448979593,1,0.7904122778426745,0.8345316311433326,nan,nan,0.11302378885599486 +204.08163265306123,1,0.7743165421613408,0.8511262472797302,nan,nan,0.12318805546636806 +255.10204081632654,1,0.753249548414535,0.8630503873758559,nan,nan,0.1361800248346027 +306.12244897959187,1,0.7268712100315663,0.8699046823746542,nan,nan,0.1530226526335016 +357.14285714285717,1,0.694762808612201,0.871183254587578,nan,nan,0.17526767302647034 +408.16326530612247,1,0.6564196199390899,0.8662543021471459,nan,nan,0.2053375622554458 +459.18367346938777,1,0.6112440678303533,0.8543369207025031,nan,nan,0.2471342724034155 +510.2040816326531,1,0.5585412933875817,0.8344743818340995,nan,nan,0.30714302736709537 +561.2244897959184,1,0.49752057854752274,0.8055051180185213,nan,nan,0.39642605066463393 +612.2448979591837,1,0.42730869340142713,0.7660347500990113,nan,nan,0.5339941531274536 +663.265306122449,1,0.3469855704957162,0.7144164850051454,nan,nan,0.7510921949147653 +714.2857142857143,1,0.25565954973242283,0.6487544841765712,nan,nan,1.089659065489047 +765.3061224489796,1,0.15260947588849744,0.5669573266103906,nan,nan,1.5668452007999316 +816.3265306122449,1,0.037533685832437026,0.4668885487044668,nan,nan,2.062626571880702 +867.3469387755102,1,-0.08904312701616794,0.34668838046440426,nan,nan,2.2566005333932724 +918.3673469387755,1,-0.22515974344986556,0.2053664293775624,nan,nan,1.9742762557365308 +969.3877551020408,1,-0.36661560366594714,0.04375763921362707,nan,nan,1.4724201265196528 +1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,nan,nan,1.0378267478370446 +1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,nan,nan,0.7457820123781923 +1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,nan,nan,0.5692403465256193 +1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,nan,nan,0.3677023181237994 +1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,nan,nan,0.310805460038203 +1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,nan,nan,0.41184810403842464 +1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,nan,nan,0.4267037403343691 +1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,nan,nan,0.45692118233449847 +1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,nan,nan,0.4835731029852539 +1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,nan,nan,0.47696040679432183 +1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,nan,nan,0.4173135946760028 +1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,nan,nan,0.32665907139175554 +1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,nan,nan,0.24572306795533427 +1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,nan,nan,0.189864891641152 +1734.6938775510205,1,0.246302320136533,-0.18222582266793214,nan,nan,0.1541293228877271 +1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,nan,nan,0.1309447698576423 +1836.734693877551,1,0.4082247422781462,0.007108706138521259,nan,nan,0.11532666162626576 +1887.7551020408164,1,0.47580638683140664,0.09358245909783232,nan,nan,0.10445257538616691 +1938.7755102040817,1,0.5354357374912007,0.17456889152858981,nan,nan,0.09672357063471888 +1989.795918367347,1,0.5877970769939762,0.25022587899384413,nan,nan,0.09119647427470219 +2040.8163265306123,1,0.6335080358663165,0.3207715550593558,nan,nan,0.08729085845318117 +2091.8367346938776,1,0.673110446052766,0.3864454566166559,nan,nan,0.08463812457943871 +2142.857142857143,1,0.707069071798058,0.44748437788656903,nan,nan,0.0830008016213394 +2193.877551020408,1,0.7357742511645264,0.5041075168486121,nan,nan,0.08222762558027762 +2244.8979591836737,1,0.75954615418438,0.5565073629841016,nan,nan,0.08222788911616931 +2295.918367346939,1,0.7786393591852832,0.6048440073048987,nan,nan,0.08295689503656585 +2346.9387755102043,1,0.7932470300017178,0.6492413535736495,nan,nan,0.0844083805697474 +2397.9591836734694,1,0.8035043081850315,0.6897842167972191,nan,nan,0.08661187268900492 +2448.979591836735,1,0.8094907181927744,0.7265156142343624,nan,nan,0.08963415615760459 +2500.0,1,0.811231477773281,0.7594337387611084,nan,nan,0.09358490399640412 From ddd6820a3b37813d7eaa7292bb13c51536dedfbe Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 11 Oct 2024 13:50:27 -0500 Subject: [PATCH 31/44] changed errors to 0.01 --- .../example_data/orbital_data_with_id.csv | 100 +++++++++--------- 1 file changed, 50 insertions(+), 50 deletions(-) diff --git a/orbitize/example_data/orbital_data_with_id.csv b/orbitize/example_data/orbital_data_with_id.csv index 5560a0f3..2c8c36e8 100644 --- a/orbitize/example_data/orbital_data_with_id.csv +++ b/orbitize/example_data/orbital_data_with_id.csv @@ -1,51 +1,51 @@ epoch,object,raoff,decoff,raoff_err,decoff_err,brightness -0.0,1,0.811231477773281,0.7594337387611091,nan,nan,0.09358490399640425 -51.02040816326531,1,0.8086976534668815,0.7884882174763149,nan,nan,0.09862733963920862 -102.04081632653062,1,0.8018051182543783,0.8135752998448389,nan,nan,0.10499582237388896 -153.06122448979593,1,0.7904122778426745,0.8345316311433326,nan,nan,0.11302378885599486 -204.08163265306123,1,0.7743165421613408,0.8511262472797302,nan,nan,0.12318805546636806 -255.10204081632654,1,0.753249548414535,0.8630503873758559,nan,nan,0.1361800248346027 -306.12244897959187,1,0.7268712100315663,0.8699046823746542,nan,nan,0.1530226526335016 -357.14285714285717,1,0.694762808612201,0.871183254587578,nan,nan,0.17526767302647034 -408.16326530612247,1,0.6564196199390899,0.8662543021471459,nan,nan,0.2053375622554458 -459.18367346938777,1,0.6112440678303533,0.8543369207025031,nan,nan,0.2471342724034155 -510.2040816326531,1,0.5585412933875817,0.8344743818340995,nan,nan,0.30714302736709537 -561.2244897959184,1,0.49752057854752274,0.8055051180185213,nan,nan,0.39642605066463393 -612.2448979591837,1,0.42730869340142713,0.7660347500990113,nan,nan,0.5339941531274536 -663.265306122449,1,0.3469855704957162,0.7144164850051454,nan,nan,0.7510921949147653 -714.2857142857143,1,0.25565954973242283,0.6487544841765712,nan,nan,1.089659065489047 -765.3061224489796,1,0.15260947588849744,0.5669573266103906,nan,nan,1.5668452007999316 -816.3265306122449,1,0.037533685832437026,0.4668885487044668,nan,nan,2.062626571880702 -867.3469387755102,1,-0.08904312701616794,0.34668838046440426,nan,nan,2.2566005333932724 -918.3673469387755,1,-0.22515974344986556,0.2053664293775624,nan,nan,1.9742762557365308 -969.3877551020408,1,-0.36661560366594714,0.04375763921362707,nan,nan,1.4724201265196528 -1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,nan,nan,1.0378267478370446 -1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,nan,nan,0.7457820123781923 -1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,nan,nan,0.5692403465256193 -1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,nan,nan,0.3677023181237994 -1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,nan,nan,0.310805460038203 -1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,nan,nan,0.41184810403842464 -1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,nan,nan,0.4267037403343691 -1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,nan,nan,0.45692118233449847 -1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,nan,nan,0.4835731029852539 -1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,nan,nan,0.47696040679432183 -1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,nan,nan,0.4173135946760028 -1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,nan,nan,0.32665907139175554 -1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,nan,nan,0.24572306795533427 -1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,nan,nan,0.189864891641152 -1734.6938775510205,1,0.246302320136533,-0.18222582266793214,nan,nan,0.1541293228877271 -1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,nan,nan,0.1309447698576423 -1836.734693877551,1,0.4082247422781462,0.007108706138521259,nan,nan,0.11532666162626576 -1887.7551020408164,1,0.47580638683140664,0.09358245909783232,nan,nan,0.10445257538616691 -1938.7755102040817,1,0.5354357374912007,0.17456889152858981,nan,nan,0.09672357063471888 -1989.795918367347,1,0.5877970769939762,0.25022587899384413,nan,nan,0.09119647427470219 -2040.8163265306123,1,0.6335080358663165,0.3207715550593558,nan,nan,0.08729085845318117 -2091.8367346938776,1,0.673110446052766,0.3864454566166559,nan,nan,0.08463812457943871 -2142.857142857143,1,0.707069071798058,0.44748437788656903,nan,nan,0.0830008016213394 -2193.877551020408,1,0.7357742511645264,0.5041075168486121,nan,nan,0.08222762558027762 -2244.8979591836737,1,0.75954615418438,0.5565073629841016,nan,nan,0.08222788911616931 -2295.918367346939,1,0.7786393591852832,0.6048440073048987,nan,nan,0.08295689503656585 -2346.9387755102043,1,0.7932470300017178,0.6492413535736495,nan,nan,0.0844083805697474 -2397.9591836734694,1,0.8035043081850315,0.6897842167972191,nan,nan,0.08661187268900492 -2448.979591836735,1,0.8094907181927744,0.7265156142343624,nan,nan,0.08963415615760459 -2500.0,1,0.811231477773281,0.7594337387611084,nan,nan,0.09358490399640412 +0.0,1,0.811231477773281,0.7594337387611091,0.01,0.01,0.09358490399640425 +51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.01,0.01,0.09862733963920862 +102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.01,0.01,0.10499582237388896 +153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.01,0.01,0.11302378885599486 +204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.01,0.01,0.12318805546636806 +255.10204081632654,1,0.753249548414535,0.8630503873758559,0.01,0.01,0.1361800248346027 +306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.01,0.01,0.1530226526335016 +357.14285714285717,1,0.694762808612201,0.871183254587578,0.01,0.01,0.17526767302647034 +408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.01,0.01,0.2053375622554458 +459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.01,0.01,0.2471342724034155 +510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.01,0.01,0.30714302736709537 +561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.01,0.01,0.39642605066463393 +612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.01,0.01,0.5339941531274536 +663.265306122449,1,0.3469855704957162,0.7144164850051454,0.01,0.01,0.7510921949147653 +714.2857142857143,1,0.25565954973242283,0.6487544841765712,0.01,0.01,1.089659065489047 +765.3061224489796,1,0.15260947588849744,0.5669573266103906,0.01,0.01,1.5668452007999316 +816.3265306122449,1,0.037533685832437026,0.4668885487044668,0.01,0.01,2.062626571880702 +867.3469387755102,1,-0.08904312701616794,0.34668838046440426,0.01,0.01,2.2566005333932724 +918.3673469387755,1,-0.22515974344986556,0.2053664293775624,0.01,0.01,1.9742762557365308 +969.3877551020408,1,-0.36661560366594714,0.04375763921362707,0.01,0.01,1.4724201265196528 +1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,0.01,0.01,1.0378267478370446 +1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.01,0.01,0.7457820123781923 +1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.01,0.01,0.5692403465256193 +1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.01,0.01,0.3677023181237994 +1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.01,0.01,0.310805460038203 +1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.01,0.01,0.41184810403842464 +1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.01,0.01,0.4267037403343691 +1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.01,0.01,0.45692118233449847 +1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.01,0.01,0.4835731029852539 +1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.01,0.01,0.47696040679432183 +1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.01,0.01,0.4173135946760028 +1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.01,0.01,0.32665907139175554 +1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.01,0.01,0.24572306795533427 +1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.01,0.01,0.189864891641152 +1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.01,0.01,0.1541293228877271 +1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.01,0.01,0.1309447698576423 +1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.01,0.01,0.11532666162626576 +1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.01,0.01,0.10445257538616691 +1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.01,0.01,0.09672357063471888 +1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.01,0.01,0.09119647427470219 +2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.01,0.01,0.08729085845318117 +2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.01,0.01,0.08463812457943871 +2142.857142857143,1,0.707069071798058,0.44748437788656903,0.01,0.01,0.0830008016213394 +2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.01,0.01,0.08222762558027762 +2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.01,0.01,0.08222788911616931 +2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.01,0.01,0.08295689503656585 +2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.01,0.01,0.0844083805697474 +2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.01,0.01,0.08661187268900492 +2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.01,0.01,0.08963415615760459 +2500.0,1,0.811231477773281,0.7594337387611084,0.01,0.01,0.09358490399640412 From 4d14d011905a6f4c5db756b6bc1a89f94bd210ce Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 18 Oct 2024 13:46:48 -0500 Subject: [PATCH 32/44] mcmc test --- mcmc_test.py | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 mcmc_test.py diff --git a/mcmc_test.py b/mcmc_test.py new file mode 100644 index 00000000..fbafc358 --- /dev/null +++ b/mcmc_test.py @@ -0,0 +1,40 @@ +import numpy as np + +import orbitize +from orbitize import driver +import multiprocessing as mp + +filename = "{}/GJ504.csv".format(orbitize.DATADIR) + +# system parameters +num_secondary_bodies = 1 +total_mass = 1.75 # [Msol] +plx = 51.44 # [mas] +mass_err = 0.05 # [Msol] +plx_err = 0.12 # [mas] + +# MCMC parameters +num_temps = 5 +num_walkers = 20 +num_threads = 2 # or a different number if you prefer, mp.cpu_count() for example + + +my_driver = driver.Driver( + filename, + "MCMC", + num_secondary_bodies, + total_mass, + plx, + mass_err=mass_err, + plx_err=plx_err, + mcmc_kwargs={ + "num_temps": num_temps, + "num_walkers": num_walkers, + "num_threads": num_threads, + }, +) +total_orbits = 6000 # number of steps x number of walkers (at lowest temperature) +burn_steps = 10 # steps to burn in per walker +thin = 2 # only save every 2nd step + +my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) \ No newline at end of file From f7c2208131557cabc97fa922c43a7d521bdcb445 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Thu, 31 Oct 2024 11:52:57 -0500 Subject: [PATCH 33/44] trying to run orbitize plot just with simulated ra and dec data --- mcmc_test.py | 21 ++++++-- .../example_data/simulated_ra_dec_data.csv | 51 +++++++++++++++++++ 2 files changed, 67 insertions(+), 5 deletions(-) create mode 100644 orbitize/example_data/simulated_ra_dec_data.csv diff --git a/mcmc_test.py b/mcmc_test.py index fbafc358..46f0ae78 100644 --- a/mcmc_test.py +++ b/mcmc_test.py @@ -4,7 +4,7 @@ from orbitize import driver import multiprocessing as mp -filename = "{}/GJ504.csv".format(orbitize.DATADIR) +filename = "{}/simulated_ra_dec_data.csv".format(orbitize.DATADIR) # system parameters num_secondary_bodies = 1 @@ -33,8 +33,19 @@ "num_threads": num_threads, }, ) -total_orbits = 6000 # number of steps x number of walkers (at lowest temperature) -burn_steps = 10 # steps to burn in per walker -thin = 2 # only save every 2nd step +if __name__ == '__main__': + + total_orbits = 6000 # number of steps x number of walkers (at lowest temperature) + burn_steps = 10 # steps to burn in per walker + thin = 2 # only save every 2nd step + + my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + + +if my_driver.sampler.results is not None: + print("Sampler results are available.") + print("Number of orbits:", len(my_driver.sampler.results.post)) +else: + print("No sampler results available. Sampler may not have run correctly.") + -my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) \ No newline at end of file diff --git a/orbitize/example_data/simulated_ra_dec_data.csv b/orbitize/example_data/simulated_ra_dec_data.csv new file mode 100644 index 00000000..d852e0c5 --- /dev/null +++ b/orbitize/example_data/simulated_ra_dec_data.csv @@ -0,0 +1,51 @@ +epoch,object,raoff,decoff,raoff_err,decoff_err +0.0,1,0.811231477773281,0.7594337387611091,0.01,0.01 +51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.01,0.01 +102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.01,0.01 +153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.01,0.01 +204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.01,0.01 +255.10204081632654,1,0.753249548414535,0.8630503873758559,0.01,0.01 +306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.01,0.01 +357.14285714285717,1,0.694762808612201,0.871183254587578,0.01,0.01 +408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.01,0.01 +459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.01,0.01 +510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.01,0.01 +561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.01,0.01 +612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.01,0.01 +663.265306122449,1,0.3469855704957162,0.7144164850051454,0.01,0.01 +714.2857142857143,1,0.25565954973242283,0.6487544841765712,0.01,0.01 +765.3061224489796,1,0.15260947588849744,0.5669573266103906,0.01,0.01 +816.3265306122449,1,0.037533685832437026,0.4668885487044668,0.01,0.01 +867.3469387755102,1,-0.08904312701616794,0.34668838046440426,0.01,0.01 +918.3673469387755,1,-0.22515974344986556,0.2053664293775624,0.01,0.01 +969.3877551020408,1,-0.36661560366594714,0.04375763921362707,0.01,0.01 +1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,0.01,0.01 +1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.01,0.01 +1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.01,0.01 +1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.01,0.01 +1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.01,0.01 +1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.01,0.01 +1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.01,0.01 +1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.01,0.01 +1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.01,0.01 +1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.01,0.01 +1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.01,0.01 +1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.01,0.01 +1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.01,0.01 +1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.01,0.01 +1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.01,0.01 +1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.01,0.01 +1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.01,0.01 +1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.01,0.01 +1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.01,0.01 +1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.01,0.01 +2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.01,0.01 +2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.01,0.01 +2142.857142857143,1,0.707069071798058,0.44748437788656903,0.01,0.01 +2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.01,0.01 +2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.01,0.01 +2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.01,0.01 +2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.01,0.01 +2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.01,0.01 +2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.01,0.01 +2500.0,1,0.811231477773281,0.7594337387611084,0.01,0.01 From 8890f96f3528978f25e401ff3be6b2f1dabd0424 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Thu, 31 Oct 2024 12:11:20 -0500 Subject: [PATCH 34/44] Trying to find what's wrong in simulated ra and dec data file --- mcmc_test.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/mcmc_test.py b/mcmc_test.py index 46f0ae78..378826c2 100644 --- a/mcmc_test.py +++ b/mcmc_test.py @@ -42,10 +42,22 @@ my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) +try: + my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + print("Sampler ran successfully.") +except Exception as e: + print(f"Error during sampling: {e}") + if my_driver.sampler.results is not None: print("Sampler results are available.") - print("Number of orbits:", len(my_driver.sampler.results.post)) + print("Results attributes:", dir(my_driver.sampler.results)) + + if my_driver.sampler.results.post is not None: + print("Posterior samples found.") + else: + print("No posteriors generated; `post` is None.") else: - print("No sampler results available. Sampler may not have run correctly.") + print("No sampler results available.") + From e70e82a931ae188746da93d6d6915f8ec621e0f0 Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Mon, 3 Feb 2025 21:06:18 -0600 Subject: [PATCH 35/44] changed sep pa end year --- mcmc_test.py | 59 ++++++++++++++++++++++++++++------------------ orbitize/kepler.py | 2 +- orbitize/plot.py | 4 ++-- 3 files changed, 39 insertions(+), 26 deletions(-) diff --git a/mcmc_test.py b/mcmc_test.py index 378826c2..30ba712d 100644 --- a/mcmc_test.py +++ b/mcmc_test.py @@ -14,9 +14,9 @@ plx_err = 0.12 # [mas] # MCMC parameters -num_temps = 5 -num_walkers = 20 -num_threads = 2 # or a different number if you prefer, mp.cpu_count() for example +num_temps = 20 +num_walkers = 1000 +num_threads = 20 # or a different number if you prefer, mp.cpu_count() for example my_driver = driver.Driver( @@ -31,33 +31,46 @@ "num_temps": num_temps, "num_walkers": num_walkers, "num_threads": num_threads, - }, + }, system_kwargs={"restrict_angle_ranges": True}, ) + + if __name__ == '__main__': - total_orbits = 6000 # number of steps x number of walkers (at lowest temperature) - burn_steps = 10 # steps to burn in per walker - thin = 2 # only save every 2nd step + total_orbits = 50000000 # number of steps x number of walkers (at lowest temperature) + burn_steps = 10000 # steps to burn in per walker + thin = 10 # only save every 2nd step my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + epochs = my_driver.system.data_table["epoch"] -try: - my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) - print("Sampler ran successfully.") -except Exception as e: - print(f"Error during sampling: {e}") - -if my_driver.sampler.results is not None: - print("Sampler results are available.") - print("Results attributes:", dir(my_driver.sampler.results)) - - if my_driver.sampler.results.post is not None: - print("Posterior samples found.") - else: - print("No posteriors generated; `post` is None.") -else: - print("No sampler results available.") + orbit_plot_fig = my_driver.sampler.results.plot_orbits( + object_to_plot=1, # Plot orbits for the first (and only, in this case) companion + num_orbits_to_plot=100, # Will plot 100 randomly selected orbits of this companion + start_mjd=epochs[0], # Minimum MJD for colorbar (here we choose first data epoch) + sep_pa_end_year=1870.0 + ) + orbit_plot_fig.savefig( + "my_orbit_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + + sma_chains, ecc_chains = my_driver.sampler.examine_chains( + param_list=["sma1", "ecc1"], n_walkers=5 + ) + + from orbitize import results + + + hdf5_filename = "my_posterior.hdf5" + import os +# To avoid weird behaviours, delete saved file if it already exists from a previous run of this notebook + if os.path.isfile(hdf5_filename): + os.remove(hdf5_filename) + my_driver.sampler.results.save_results(hdf5_filename) + loaded_results = results.Results() # Create blank results object for loading + loaded_results.load_results(hdf5_filename) + diff --git a/orbitize/kepler.py b/orbitize/kepler.py index c3bae2be..60e4d5af 100644 --- a/orbitize/kepler.py +++ b/orbitize/kepler.py @@ -58,7 +58,7 @@ def times2trueanom_and_eccanom( use_gpu=False, ): - print("hi Farrah!!") + n_orbs = np.size(sma) # num sets of input orbital parameters n_dates = np.size(epochs) # number of dates to compute offsets and vz diff --git a/orbitize/plot.py b/orbitize/plot.py index a02cfc1e..1177115f 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -173,7 +173,7 @@ def plot_orbits( show_colorbar=True, cmap=cmap, sep_pa_color="lightgrey", - sep_pa_end_year=2025.0, + sep_pa_end_year=1875.0, cbar_param="Epoch [year]", mod180=False, rv_time_series=False, @@ -1077,7 +1077,7 @@ def plot_residuals( num_orbits_to_plot=100, num_epochs_to_plot=100, sep_pa_color="lightgrey", - sep_pa_end_year=2025.0, + sep_pa_end_year=1925.0, cbar_param="Epoch [year]", mod180=False, ): From ffcaacefa7ed54846295c947be7f86187a508f7d Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Fri, 14 Feb 2025 13:53:14 -0800 Subject: [PATCH 36/44] brightness calc incorp in mcmc --- .../example_data/orbital_data_with_id.csv | 102 +++++++++--------- .../example_data/reflected_light_example.csv | 20 ++-- orbitize/read_input.py | 37 +++---- orbitize/system.py | 7 +- tests/test_brightness.py | 59 +++++----- 5 files changed, 107 insertions(+), 118 deletions(-) diff --git a/orbitize/example_data/orbital_data_with_id.csv b/orbitize/example_data/orbital_data_with_id.csv index 2c8c36e8..c8f15a9d 100644 --- a/orbitize/example_data/orbital_data_with_id.csv +++ b/orbitize/example_data/orbital_data_with_id.csv @@ -1,51 +1,51 @@ -epoch,object,raoff,decoff,raoff_err,decoff_err,brightness -0.0,1,0.811231477773281,0.7594337387611091,0.01,0.01,0.09358490399640425 -51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.01,0.01,0.09862733963920862 -102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.01,0.01,0.10499582237388896 -153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.01,0.01,0.11302378885599486 -204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.01,0.01,0.12318805546636806 -255.10204081632654,1,0.753249548414535,0.8630503873758559,0.01,0.01,0.1361800248346027 -306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.01,0.01,0.1530226526335016 -357.14285714285717,1,0.694762808612201,0.871183254587578,0.01,0.01,0.17526767302647034 -408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.01,0.01,0.2053375622554458 -459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.01,0.01,0.2471342724034155 -510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.01,0.01,0.30714302736709537 -561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.01,0.01,0.39642605066463393 -612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.01,0.01,0.5339941531274536 -663.265306122449,1,0.3469855704957162,0.7144164850051454,0.01,0.01,0.7510921949147653 -714.2857142857143,1,0.25565954973242283,0.6487544841765712,0.01,0.01,1.089659065489047 -765.3061224489796,1,0.15260947588849744,0.5669573266103906,0.01,0.01,1.5668452007999316 -816.3265306122449,1,0.037533685832437026,0.4668885487044668,0.01,0.01,2.062626571880702 -867.3469387755102,1,-0.08904312701616794,0.34668838046440426,0.01,0.01,2.2566005333932724 -918.3673469387755,1,-0.22515974344986556,0.2053664293775624,0.01,0.01,1.9742762557365308 -969.3877551020408,1,-0.36661560366594714,0.04375763921362707,0.01,0.01,1.4724201265196528 -1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,0.01,0.01,1.0378267478370446 -1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.01,0.01,0.7457820123781923 -1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.01,0.01,0.5692403465256193 -1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.01,0.01,0.3677023181237994 -1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.01,0.01,0.310805460038203 -1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.01,0.01,0.41184810403842464 -1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.01,0.01,0.4267037403343691 -1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.01,0.01,0.45692118233449847 -1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.01,0.01,0.4835731029852539 -1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.01,0.01,0.47696040679432183 -1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.01,0.01,0.4173135946760028 -1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.01,0.01,0.32665907139175554 -1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.01,0.01,0.24572306795533427 -1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.01,0.01,0.189864891641152 -1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.01,0.01,0.1541293228877271 -1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.01,0.01,0.1309447698576423 -1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.01,0.01,0.11532666162626576 -1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.01,0.01,0.10445257538616691 -1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.01,0.01,0.09672357063471888 -1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.01,0.01,0.09119647427470219 -2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.01,0.01,0.08729085845318117 -2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.01,0.01,0.08463812457943871 -2142.857142857143,1,0.707069071798058,0.44748437788656903,0.01,0.01,0.0830008016213394 -2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.01,0.01,0.08222762558027762 -2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.01,0.01,0.08222788911616931 -2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.01,0.01,0.08295689503656585 -2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.01,0.01,0.0844083805697474 -2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.01,0.01,0.08661187268900492 -2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.01,0.01,0.08963415615760459 -2500.0,1,0.811231477773281,0.7594337387611084,0.01,0.01,0.09358490399640412 +epoch,object,raoff,decoff,raoff_err,decoff_err,brightness,brightness_err +0.0,1,0.811231477773281,0.7594337387611091,0.01,0.01,0.09358490399640425,0.01 +51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.01,0.01,0.09862733963920862,0.01 +102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.01,0.01,0.10499582237388896,0.01 +153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.01,0.01,0.11302378885599486,0.01 +204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.01,0.01,0.12318805546636806,0.01 +255.10204081632654,1,0.753249548414535,0.8630503873758559,0.01,0.01,0.1361800248346027,0.01 +306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.01,0.01,0.1530226526335016,0.01 +357.14285714285717,1,0.694762808612201,0.871183254587578,0.01,0.01,0.17526767302647034,0.01 +408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.01,0.01,0.2053375622554458,0.01 +459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.01,0.01,0.2471342724034155,0.01 +510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.01,0.01,0.30714302736709537,0.01 +561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.01,0.01,0.39642605066463393,0.01 +612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.01,0.01,0.5339941531274536,0.01 +663.265306122449,1,0.3469855704957162,0.7144164850051454,0.01,0.01,0.7510921949147653,0.01 +714.2857142857143,1,0.25565954973242283,0.6487544841765712,0.01,0.01,1.089659065489047,0.01 +765.3061224489796,1,0.15260947588849744,0.5669573266103906,0.01,0.01,1.5668452007999316,0.01 +816.3265306122449,1,0.037533685832437026,0.4668885487044668,0.01,0.01,2.062626571880702,0.01 +867.3469387755102,1,-0.08904312701616794,0.34668838046440426,0.01,0.01,2.2566005333932724,0.01 +918.3673469387755,1,-0.22515974344986556,0.2053664293775624,0.01,0.01,1.9742762557365308,0.01 +969.3877551020408,1,-0.36661560366594714,0.04375763921362707,0.01,0.01,1.4724201265196528,0.01 +1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,0.01,0.01,1.0378267478370446,0.01 +1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.01,0.01,0.7457820123781923,0.01 +1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.01,0.01,0.5692403465256193,0.01 +1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.01,0.01,0.3677023181237994,0.01 +1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.01,0.01,0.310805460038203,0.01 +1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.01,0.01,0.41184810403842464,0.01 +1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.01,0.01,0.4267037403343691,0.01 +1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.01,0.01,0.45692118233449847,0.01 +1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.01,0.01,0.4835731029852539,0.01 +1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.01,0.01,0.47696040679432183,0.01 +1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.01,0.01,0.4173135946760028,0.01 +1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.01,0.01,0.32665907139175554,0.01 +1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.01,0.01,0.24572306795533427,0.01 +1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.01,0.01,0.189864891641152,0.01 +1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.01,0.01,0.1541293228877271,0.01 +1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.01,0.01,0.1309447698576423,0.01 +1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.01,0.01,0.11532666162626576,0.01 +1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.01,0.01,0.10445257538616691,0.01 +1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.01,0.01,0.09672357063471888,0.01 +1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.01,0.01,0.09119647427470219,0.01 +2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.01,0.01,0.08729085845318117,0.01 +2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.01,0.01,0.08463812457943871,0.01 +2142.857142857143,1,0.707069071798058,0.44748437788656903,0.01,0.01,0.0830008016213394,0.01 +2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.01,0.01,0.08222762558027762,0.01 +2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.01,0.01,0.08222788911616931,0.01 +2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.01,0.01,0.08295689503656585,0.01 +2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.01,0.01,0.0844083805697474,0.01 +2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.01,0.01,0.08661187268900492,0.01 +2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.01,0.01,0.08963415615760459,0.01 +2500.0,1,0.811231477773281,0.7594337387611084,0.01,0.01,0.09358490399640412,0.01 diff --git a/orbitize/example_data/reflected_light_example.csv b/orbitize/example_data/reflected_light_example.csv index 8f0fcec9..1d7cd32a 100644 --- a/orbitize/example_data/reflected_light_example.csv +++ b/orbitize/example_data/reflected_light_example.csv @@ -1,10 +1,10 @@ -epoch,object,raoff,decoff,raoff_err,decoff_err,brightness -57298,1,253.72,92.35,2.98,2.85,0.5 -57606,1,236.63,127.94,9.77,9.18,0.5 -57645,1,234.52,123.39,1.79,1.03,0.5 -57946,1,210.76,152.09,1.94,1.88,0.5 -58276,1,167.49,180.87,1.61,16.97, -58287,1,177.67,174.6,1.67,1.67,0.5 -58365,1,165.7,185.33,3.28,3.66,0.5 -58368,1,170.38,185.94,2.52,2.74,0.5 -58414,1,161.64,176.21,13.6,14.31, \ No newline at end of file +epoch,object,raoff,decoff,raoff_err,decoff_err,brightness,brightness_err +57298,1,253.72,92.35,2.98,2.85,0.5,0.1 +57606,1,236.63,127.94,9.77,9.18,0.5,0.1 +57645,1,234.52,123.39,1.79,1.03,0.5,0.1 +57946,1,210.76,152.09,1.94,1.88,0.5,0.1 +58276,1,167.49,180.87,1.61,16.97,, +58287,1,177.67,174.6,1.67,1.67,0.5,0.1 +58365,1,165.7,185.33,3.28,3.66,0.5,0.1 +58368,1,170.38,185.94,2.52,2.74,0.5,0.1 +58414,1,161.64,176.21,13.6,14.31,, \ No newline at end of file diff --git a/orbitize/read_input.py b/orbitize/read_input.py index 9927c706..cd6f6a9a 100644 --- a/orbitize/read_input.py +++ b/orbitize/read_input.py @@ -103,9 +103,8 @@ def read_file(filename): "quant12_corr", "quant_type", "instrument", - "brightness", ), - dtype=(float, int, float, float, float, float, float, "S5", "S5", float), + dtype=(float, int, float, float, float, float, float, "S5", "S5"), ) # read file @@ -251,12 +250,6 @@ def read_file(filename): if "quant12_corr" not in input_table.keys(): default_corrs = np.nan * np.ones(len(input_table)) input_table.add_column(default_corrs, name="quant12_corr") - if "brightness" not in input_table.keys(): - default_brightness = np.nan * np.ones(len(input_table)) - input_table.add_column(default_brightness, name="brightness") - have_brightness = np.zeros(num_measurements, dtype=bool) - else: - have_brightness = np.ones(num_measurements, dtype=bool) if "instrument" not in input_table.keys(): default_insts = [] for this_quant_type in input_table["quant_type"]: @@ -289,18 +282,14 @@ def read_file(filename): MJD = row["epoch"] - 2400000.5 else: MJD = row["epoch"] - if have_brightness[index]: - brightness = row["brightness"] - else: - brightness = np.nan # check that "object" is an integer (instead of ABC/bcd) if not isinstance(row["object"], (int, np.int32, np.int64)): raise Exception("Invalid object ID. Object IDs must be integers.") - # determine input quantity type (RA/DEC, SEP/PA, or RV) + # determine input quantity type (RA/DEC, SEP/PA, RV, or BRIGHTNESS) if orbitize_style: - if row["quant_type"] == "rv": # special format for rv rows + if row["quant_type"] == "rv" or row["quant_type"] == 'brightness': # special format for rv rows output_table.add_row( [ MJD, @@ -312,7 +301,6 @@ def read_file(filename): None, row["quant_type"], row["instrument"], - None, ] ) @@ -341,7 +329,6 @@ def read_file(filename): quant12_corr, row["quant_type"], row["instrument"], - brightness, ] ) else: # catch wrong formats @@ -383,7 +370,6 @@ def read_file(filename): quant12_corr, "radec", this_inst, - brightness, ] ) @@ -417,7 +403,6 @@ def read_file(filename): quant12_corr, "seppa", this_inst, - brightness, ] ) @@ -434,7 +419,6 @@ def read_file(filename): None, "rv", row["instrument"], - brightness, ] ) else: @@ -450,9 +434,22 @@ def read_file(filename): None, "rv", "defrv", - brightness, ] ) + if have_brightness[index]: + output_table.add_row( + [ + MJD, + row["object"], + row["brightness"], + row["brightness_err"], + None, + None, + None, + "brightness", + "defbr", + ] + ) return output_table diff --git a/orbitize/system.py b/orbitize/system.py index d2ed093c..55e8a97a 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -105,7 +105,7 @@ def __init__( radec_indices = np.where(self.data_table["quant_type"] == "radec") seppa_indices = np.where(self.data_table["quant_type"] == "seppa") - brightness_indices = np.where(~np.isnan(self.data_table["brightness"])) + brightness_indices = np.where(self.data_table["quant_type"] == "brightness") if len(radec_indices[0]) == 0 and len(seppa_indices[0]) == 0: self.fit_astrometry = False @@ -387,6 +387,7 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): dec_perturb = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) vz = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) + brightness_out = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) # mass/mtot used to compute each Keplerian orbit will be needed later to compute perturbations if self.track_planet_perturbs: @@ -510,11 +511,13 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): raoff = np.array([raoff]) decoff = np.array([decoff]) vz_i = np.array([vz_i]) + brightness_out = np.array([brightness]) # add Keplerian ra/deoff for this body to storage arrays ra_kepler[:, body_num, :] = np.reshape(raoff, (n_epochs, n_orbits)) dec_kepler[:, body_num, :] = np.reshape(decoff, (n_epochs, n_orbits)) vz[:, body_num, :] = np.reshape(vz_i, (n_epochs, n_orbits)) + brightness_out[:, body_num, :] = np.reshape(brightness, (n_epochs, n_orbits)) # vz_i is the ith companion radial velocity if self.fit_secondary_mass: @@ -586,7 +589,7 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): else: return raoff, deoff, vz else: - return raoff, deoff, vz, brightness + return raoff, deoff, vz, brightness_out def compute_model(self, params_arr, use_rebound=False): """ diff --git a/tests/test_brightness.py b/tests/test_brightness.py index 60091cd1..7a3f3d94 100644 --- a/tests/test_brightness.py +++ b/tests/test_brightness.py @@ -2,17 +2,16 @@ Tests for the reflected-light calculation in system.compute_all_orbits """ -from orbitize import system, sampler +from orbitize import system, sampler, plot from orbitize import DATADIR, read_input import os import numpy as np import pandas as pd import matplotlib.pyplot as plt +from astropy.time import Time -print('hello') - -def test_brightness_calculation(): +def test_brightness_calculation(make_plot=False): num_secondary_bodies = 1 @@ -27,8 +26,6 @@ def test_brightness_calculation(): test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) - print(test_system.param_idx) - params = np.array( [ 10.0, @@ -44,41 +41,29 @@ def test_brightness_calculation(): ra, dec, vz, brightness = test_system.compute_all_orbits(params) - plt.figure() - plt.scatter(times, brightness) - plt.xlabel("Time [dy]", fontsize=18) - plt.ylabel("Brightness", fontsize=18) - plt.savefig("Test_brightness2.png") + if make_plot: + plt.figure() + plt.scatter(times, brightness) + plt.xlabel("Time [dy]", fontsize=18) + plt.ylabel("Brightness", fontsize=18) + plt.savefig("Test_brightness2.png") def test_read_input_with_brightness(): num_secondary_bodies = 1 - input_file = os.path.join(DATADIR, "test_orbital_data.csv") + input_file = os.path.join(DATADIR, "reflected_light_example.csv") data_table = read_input.read_file(input_file) times = data_table["epoch"].value - brightness_values = data_table["brightness"].value print(data_table) # TODO (Farrah): add a test that asserts the brightness column of the data table is # what you expect (hint: check in the reflected_light_example.csv to see what # the brightness values should be -def test_assert_brightness(): - - num_secondary_bodies = 1 - - input_file = os.path.join(DATADIR, "reflected_light_example.csv") - - data_table = read_input.read_file(input_file) - - brightness_values = data_table["brightness"].value - assert 'brightness' in data_table.columns, "Brightness column does not exist" - #print("The assert works!") - #print (brightness_values) def test_assert_nan(): @@ -113,30 +98,34 @@ def test_compute_posteriors(): np.radians(90), # pan 0.0, # tau 51.5, # plx - 1.75, # stellar maxx + 1.75, # stellar mass ] ) epochs = np.linspace(0, 365*30, int(1e3)) ra, dec, vz, brightness = test_system.compute_all_orbits(params_arr, epochs=epochs) - fig, ax = plt.subplots(2, 1, figsize=(5,10)) + # fig, ax = plt.subplots(2, 1, figsize=(5,10)) - ax[0].scatter(epochs, brightness, color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) - ax[1].scatter(ra[:,1,:], dec[:,1,:], color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) + # ax[0].scatter(epochs, brightness, color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) + # ax[1].scatter(ra[:,1,:], dec[:,1,:], color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) - ax[1].axis('equal') - plt.savefig('visual4farrah.png') + # ax[1].axis('equal') + # plt.savefig('visual4farrah.png') model = test_system.compute_model(params_arr) - # print(model) - test_mcmc = sampler.MCMC(test_system, 1, 50, num_threads=1) + test_mcmc = sampler.MCMC(test_system, 1, 50) test_mcmc.run_sampler(10) + # fig = plot.plot_orbits(test_mcmc.results, start_mjd=Time(2000., format='decimalyear').mjd, sep_pa_end_year=2500.) + + # plt.savefig('foo.png') + + + if __name__ == "__main__": #test_brightness_calculation() - #test_read_input_with_brightness() - # test_assert_brightness() + # test_read_input_with_brightness() test_compute_posteriors() From e009c18386e73ca80a57d56dc5307955c8c2b0da Mon Sep 17 00:00:00 2001 From: Farrah Date: Tue, 18 Feb 2025 22:58:14 -0800 Subject: [PATCH 37/44] new test mcmc file with added brightness values --- tests/brightness_mcmc_test.py | 75 +++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 tests/brightness_mcmc_test.py diff --git a/tests/brightness_mcmc_test.py b/tests/brightness_mcmc_test.py new file mode 100644 index 00000000..e82a5d07 --- /dev/null +++ b/tests/brightness_mcmc_test.py @@ -0,0 +1,75 @@ +import numpy as np + +import orbitize +from orbitize import driver +import multiprocessing as mp + +filename = "{}/orbital_data_with_id.csv".format(orbitize.DATADIR) + +# system parameters +num_secondary_bodies = 1 +total_mass = 1.75 # [Msol] +plx = 51.44 # [mas] +mass_err = 0.05 # [Msol] +plx_err = 0.12 # [mas] + +# MCMC parameters +num_temps = 20 +num_walkers = 1000 +num_threads = 20 # or a different number if you prefer, mp.cpu_count() for example + + +my_driver = driver.Driver( + filename, + "MCMC", + num_secondary_bodies, + total_mass, + plx, + mass_err=mass_err, + plx_err=plx_err, + mcmc_kwargs={ + "num_temps": num_temps, + "num_walkers": num_walkers, + "num_threads": num_threads, + }, system_kwargs={"restrict_angle_ranges": True}, +) + + +if __name__ == '__main__': + + total_orbits = 100000000 # number of steps x number of walkers (at lowest temperature) + burn_steps = 50000 # steps to burn in per walker + thin = 10 # only save every 2nd step + + my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + + epochs = my_driver.system.data_table["epoch"] + + orbit_plot_fig = my_driver.sampler.results.plot_orbits( + object_to_plot=1, # Plot orbits for the first (and only, in this case) companion + num_orbits_to_plot=100, # Will plot 100 randomly selected orbits of this companion + start_mjd=epochs[0], # Minimum MJD for colorbar (here we choose first data epoch) + sep_pa_end_year=1870.0 + ) + orbit_plot_fig.savefig( + "my_orbit_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + + sma_chains, ecc_chains = my_driver.sampler.examine_chains( + param_list=["sma1", "ecc1"], n_walkers=5 + ) + + from orbitize import results + + + hdf5_filename = "my_posterior.hdf5" + import os + +# To avoid weird behaviours, delete saved file if it already exists from a previous run of this notebook + if os.path.isfile(hdf5_filename): + os.remove(hdf5_filename) + + my_driver.sampler.results.save_results(hdf5_filename) + loaded_results = results.Results() # Create blank results object for loading + loaded_results.load_results(hdf5_filename) + From 8e9a8e2c0df21942367b9aadcf571be5b79d0f92 Mon Sep 17 00:00:00 2001 From: Farrah Date: Mon, 24 Feb 2025 08:43:54 -0800 Subject: [PATCH 38/44] separate test for brightness mcmc --- tests/brightness_mcmc_test.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/brightness_mcmc_test.py b/tests/brightness_mcmc_test.py index e82a5d07..1e4dbeb4 100644 --- a/tests/brightness_mcmc_test.py +++ b/tests/brightness_mcmc_test.py @@ -43,6 +43,13 @@ my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + corner_plot_fig = ( + my_driver.sampler.results.plot_corner() + ) # Creates a corner plot and returns Figure object + corner_plot_fig.savefig( + "my_corner_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + epochs = my_driver.system.data_table["epoch"] orbit_plot_fig = my_driver.sampler.results.plot_orbits( From 4966b3dba0d5ca0306b567b511af6b7bf50bcfa3 Mon Sep 17 00:00:00 2001 From: Farrah Date: Tue, 4 Mar 2025 12:03:18 -0800 Subject: [PATCH 39/44] corner plot --- tests/save_chain_corner_plot.py | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 tests/save_chain_corner_plot.py diff --git a/tests/save_chain_corner_plot.py b/tests/save_chain_corner_plot.py new file mode 100644 index 00000000..9d02bd55 --- /dev/null +++ b/tests/save_chain_corner_plot.py @@ -0,0 +1,9 @@ +# Save Chains Corner Plot +import orbitize.results + +import matplotlib.pyplot as plt + +myResults = orbitize.results.Results() # create empty Results object +myResults.load_results("my_posterior.hdf5") # load from file + +corner_figure = myResults.plot_corner() \ No newline at end of file From ab6b40c6eff35118bfb63f75d1b575d5cfa93e4e Mon Sep 17 00:00:00 2001 From: Farrah Molina Date: Fri, 14 Mar 2025 14:26:34 -0500 Subject: [PATCH 40/44] new simulated ra and dec data --- HDF5_Test.py | 8 ++ .../example_data/simulated_ra_dec_data.csv | 100 +++++++++--------- 2 files changed, 58 insertions(+), 50 deletions(-) create mode 100644 HDF5_Test.py diff --git a/HDF5_Test.py b/HDF5_Test.py new file mode 100644 index 00000000..2c80b04d --- /dev/null +++ b/HDF5_Test.py @@ -0,0 +1,8 @@ +from orbitize import results +hdf5_filename = "my_posterior.hdf5" + + +loaded_results = results.Results() # Create blank results object for loading +loaded_results.load_results(hdf5_filename) + +print(loaded_results.post) \ No newline at end of file diff --git a/orbitize/example_data/simulated_ra_dec_data.csv b/orbitize/example_data/simulated_ra_dec_data.csv index d852e0c5..f70d14ff 100644 --- a/orbitize/example_data/simulated_ra_dec_data.csv +++ b/orbitize/example_data/simulated_ra_dec_data.csv @@ -1,51 +1,51 @@ epoch,object,raoff,decoff,raoff_err,decoff_err -0.0,1,0.811231477773281,0.7594337387611091,0.01,0.01 -51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.01,0.01 -102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.01,0.01 -153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.01,0.01 -204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.01,0.01 -255.10204081632654,1,0.753249548414535,0.8630503873758559,0.01,0.01 -306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.01,0.01 -357.14285714285717,1,0.694762808612201,0.871183254587578,0.01,0.01 -408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.01,0.01 -459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.01,0.01 -510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.01,0.01 -561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.01,0.01 -612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.01,0.01 -663.265306122449,1,0.3469855704957162,0.7144164850051454,0.01,0.01 -714.2857142857143,1,0.25565954973242283,0.6487544841765712,0.01,0.01 -765.3061224489796,1,0.15260947588849744,0.5669573266103906,0.01,0.01 -816.3265306122449,1,0.037533685832437026,0.4668885487044668,0.01,0.01 -867.3469387755102,1,-0.08904312701616794,0.34668838046440426,0.01,0.01 -918.3673469387755,1,-0.22515974344986556,0.2053664293775624,0.01,0.01 -969.3877551020408,1,-0.36661560366594714,0.04375763921362707,0.01,0.01 -1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,0.01,0.01 -1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.01,0.01 -1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.01,0.01 -1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.01,0.01 -1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.01,0.01 -1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.01,0.01 -1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.01,0.01 -1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.01,0.01 -1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.01,0.01 -1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.01,0.01 -1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.01,0.01 -1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.01,0.01 -1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.01,0.01 -1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.01,0.01 -1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.01,0.01 -1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.01,0.01 -1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.01,0.01 -1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.01,0.01 -1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.01,0.01 -1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.01,0.01 -2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.01,0.01 -2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.01,0.01 -2142.857142857143,1,0.707069071798058,0.44748437788656903,0.01,0.01 -2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.01,0.01 -2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.01,0.01 -2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.01,0.01 -2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.01,0.01 -2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.01,0.01 -2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.01,0.01 -2500.0,1,0.811231477773281,0.7594337387611084,0.01,0.01 +0.0,1,450.6841543184894,421.9076326450606,0.01,0.01 +51.02040816326531,1,449.2764741482675,438.04900970906385,0.01,0.01 +102.04081632653062,1,445.44728791909904,451.9862776915772,0.01,0.01 +153.06122448979593,1,439.1179321348192,463.6286839685181,0.01,0.01 +204.08163265306123,1,430.17585675630045,472.84791515540564,0.01,0.01 +255.10204081632654,1,418.4719713414083,479.4724374310311,0.01,0.01 +306.12244897959187,1,403.8173389064257,483.28037909703016,0.01,0.01 +357.14285714285717,1,385.9793381178894,483.9906969930988,0.01,0.01 +408.16326530612247,1,364.67756663282773,481.25239008174776,0.01,0.01 +459.18367346938777,1,339.5800376835296,474.63162261250176,0.01,0.01 +510.2040816326531,1,310.30071854865656,463.59687879672197,0.01,0.01 +561.2244897959184,1,276.4003214152904,447.50284334362294,0.01,0.01 +612.2448979591837,1,237.3937185563484,425.57486116611744,0.01,0.01 +663.265306122449,1,192.769761386509,396.89804722508075,0.01,0.01 +714.2857142857143,1,142.03308318467933,360.4191578758729,0.01,0.01 +765.3061224489796,1,84.78304216027635,314.9762925613281,0.01,0.01 +816.3265306122449,1,20.852047684687236,259.3825270580371,0.01,0.01 +867.3469387755102,1,-49.46840389787108,192.60465581355794,0.01,0.01 +918.3673469387755,1,-125.08874636103643,114.09246076531244,0.01,0.01 +969.3877551020408,1,-203.6753353699706,24.30979956312615,0.01,0.01 +1020.4081632653061,1,-281.1930260579934,-74.54356985421909,0.01,0.01 +1071.4285714285716,1,-351.69816254016706,-177.8056550640092,0.01,0.01 +1122.4489795918369,1,-407.813490427722,-278.09064717234304,0.01,0.01 +1173.4693877551022,1,-442.23522156712545,-366.13080690247733,0.01,0.01 +1224.4897959183675,1,-450.0052085538893,-432.9675001605461,0.01,0.01 +1275.5102040816328,1,-430.39483123995797,-472.67765040997523,0.01,0.01 +1326.530612244898,1,-387.1222955160196,-484.04061168877496,0.01,0.01 +1377.5510204081634,1,-326.77148869768655,-470.1828585947054,0.01,0.01 +1428.5714285714287,1,-256.5528105962256,-436.7513713969537,0.01,0.01 +1479.591836734694,1,-182.63510676508716,-389.9204942329679,0.01,0.01 +1530.6122448979593,1,-109.46795465403186,-335.09403393338613,0.01,0.01 +1581.6326530612246,1,-39.83805845141467,-276.40672308912804,0.01,0.01 +1632.6530612244899,1,24.75403165429632,-216.74299816205493,0.01,0.01 +1683.6734693877552,1,83.6775818280479,-157.9747940702846,0.01,0.01 +1734.6938775510205,1,136.8346222980739,-101.23656814885118,0.01,0.01 +1785.7142857142858,1,184.42422735042004,-47.16014775087199,0.01,0.01 +1836.734693877551,1,226.791523487859,3.949281188067366,0.01,0.01 +1887.7551020408164,1,264.33688157300367,51.99025505435129,0.01,0.01 +1938.7755102040817,1,297.4642986062226,96.98271751588322,0.01,0.01 +1989.795918367347,1,326.5539316633201,139.01437721880228,0.01,0.01 +2040.8163265306123,1,351.9489088146203,178.2064194774199,0.01,0.01 +2091.8367346938776,1,373.95024780709224,214.6919203425866,0.01,0.01 +2142.857142857143,1,392.8161509989211,248.60243215920502,0.01,0.01 +2193.877551020408,1,408.7634728691813,280.0597315825623,0.01,0.01 +2244.8979591836737,1,421.9700856579889,309.1707572133898,0.01,0.01 +2295.918367346939,1,432.57742176960176,336.0244485027215,0.01,0.01 +2346.9387755102043,1,440.69279444539876,360.68964087424973,0.01,0.01 +2397.9591836734694,1,446.39128232501747,383.2134537762328,0.01,0.01 +2448.979591836735,1,449.7170656626525,403.6197856857568,0.01,0.01 +2500.0,1,450.6841543184894,421.90763264506023,0.01,0.01 From a53909099593148930aa41bea031098329a8f38c Mon Sep 17 00:00:00 2001 From: Farrah Date: Sun, 16 Mar 2025 23:25:36 -0700 Subject: [PATCH 41/44] updated parallax (0.03) to match simulated ra/dec data file --- mcmc_test.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/mcmc_test.py b/mcmc_test.py index 30ba712d..2ea56883 100644 --- a/mcmc_test.py +++ b/mcmc_test.py @@ -9,7 +9,7 @@ # system parameters num_secondary_bodies = 1 total_mass = 1.75 # [Msol] -plx = 51.44 # [mas] +plx = 0.03 # [mas] mass_err = 0.05 # [Msol] plx_err = 0.12 # [mas] @@ -38,11 +38,18 @@ if __name__ == '__main__': total_orbits = 50000000 # number of steps x number of walkers (at lowest temperature) - burn_steps = 10000 # steps to burn in per walker + burn_steps = 500000 # steps to burn in per walker thin = 10 # only save every 2nd step my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + corner_plot_fig = ( + my_driver.sampler.results.plot_corner() + ) # Creates a corner plot and returns Figure object + corner_plot_fig.savefig( + "my_corner_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + epochs = my_driver.system.data_table["epoch"] orbit_plot_fig = my_driver.sampler.results.plot_orbits( @@ -72,5 +79,4 @@ my_driver.sampler.results.save_results(hdf5_filename) loaded_results = results.Results() # Create blank results object for loading loaded_results.load_results(hdf5_filename) - From 0f358a6ad500977b03ca94df58ecabf010a0434e Mon Sep 17 00:00:00 2001 From: Farrah Date: Mon, 24 Mar 2025 13:20:30 -0700 Subject: [PATCH 42/44] mcmc under tests --- mcmc_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcmc_test.py b/mcmc_test.py index 2ea56883..79c096dc 100644 --- a/mcmc_test.py +++ b/mcmc_test.py @@ -9,7 +9,7 @@ # system parameters num_secondary_bodies = 1 total_mass = 1.75 # [Msol] -plx = 0.03 # [mas] +plx = 30 # [mas] mass_err = 0.05 # [Msol] plx_err = 0.12 # [mas] From 239aff10dfa6f12f967725c1ac10cd0b2ec8ff6e Mon Sep 17 00:00:00 2001 From: Farrah Date: Fri, 6 Jun 2025 09:34:16 -0700 Subject: [PATCH 43/44] testing brightness posteriors --- HDF5_Test.py | 15 ++- .../example_data/simulated_ra_dec_data.csv | 100 +++++++++--------- 2 files changed, 63 insertions(+), 52 deletions(-) diff --git a/HDF5_Test.py b/HDF5_Test.py index 2c80b04d..4d108539 100644 --- a/HDF5_Test.py +++ b/HDF5_Test.py @@ -1,8 +1,19 @@ from orbitize import results -hdf5_filename = "my_posterior.hdf5" +import matplotlib.pyplot as plt +hdf5_filename = "my_posterior_brightness.hdf5" loaded_results = results.Results() # Create blank results object for loading loaded_results.load_results(hdf5_filename) -print(loaded_results.post) \ No newline at end of file +print(loaded_results.post) + +orbit_plot_fig = loaded_results.plot_orbits( + object_to_plot=1, + num_orbits_to_plot=100, + start_mjd=0, + sep_pa_end_year=1861 +) +# Save or show the figure +orbit_plot_fig.savefig("reloaded_short_bright_orbit_plot.png") +plt.show() \ No newline at end of file diff --git a/orbitize/example_data/simulated_ra_dec_data.csv b/orbitize/example_data/simulated_ra_dec_data.csv index f70d14ff..afc737cf 100644 --- a/orbitize/example_data/simulated_ra_dec_data.csv +++ b/orbitize/example_data/simulated_ra_dec_data.csv @@ -1,51 +1,51 @@ epoch,object,raoff,decoff,raoff_err,decoff_err -0.0,1,450.6841543184894,421.9076326450606,0.01,0.01 -51.02040816326531,1,449.2764741482675,438.04900970906385,0.01,0.01 -102.04081632653062,1,445.44728791909904,451.9862776915772,0.01,0.01 -153.06122448979593,1,439.1179321348192,463.6286839685181,0.01,0.01 -204.08163265306123,1,430.17585675630045,472.84791515540564,0.01,0.01 -255.10204081632654,1,418.4719713414083,479.4724374310311,0.01,0.01 -306.12244897959187,1,403.8173389064257,483.28037909703016,0.01,0.01 -357.14285714285717,1,385.9793381178894,483.9906969930988,0.01,0.01 -408.16326530612247,1,364.67756663282773,481.25239008174776,0.01,0.01 -459.18367346938777,1,339.5800376835296,474.63162261250176,0.01,0.01 -510.2040816326531,1,310.30071854865656,463.59687879672197,0.01,0.01 -561.2244897959184,1,276.4003214152904,447.50284334362294,0.01,0.01 -612.2448979591837,1,237.3937185563484,425.57486116611744,0.01,0.01 -663.265306122449,1,192.769761386509,396.89804722508075,0.01,0.01 -714.2857142857143,1,142.03308318467933,360.4191578758729,0.01,0.01 -765.3061224489796,1,84.78304216027635,314.9762925613281,0.01,0.01 -816.3265306122449,1,20.852047684687236,259.3825270580371,0.01,0.01 -867.3469387755102,1,-49.46840389787108,192.60465581355794,0.01,0.01 -918.3673469387755,1,-125.08874636103643,114.09246076531244,0.01,0.01 -969.3877551020408,1,-203.6753353699706,24.30979956312615,0.01,0.01 -1020.4081632653061,1,-281.1930260579934,-74.54356985421909,0.01,0.01 -1071.4285714285716,1,-351.69816254016706,-177.8056550640092,0.01,0.01 -1122.4489795918369,1,-407.813490427722,-278.09064717234304,0.01,0.01 -1173.4693877551022,1,-442.23522156712545,-366.13080690247733,0.01,0.01 -1224.4897959183675,1,-450.0052085538893,-432.9675001605461,0.01,0.01 -1275.5102040816328,1,-430.39483123995797,-472.67765040997523,0.01,0.01 -1326.530612244898,1,-387.1222955160196,-484.04061168877496,0.01,0.01 -1377.5510204081634,1,-326.77148869768655,-470.1828585947054,0.01,0.01 -1428.5714285714287,1,-256.5528105962256,-436.7513713969537,0.01,0.01 -1479.591836734694,1,-182.63510676508716,-389.9204942329679,0.01,0.01 -1530.6122448979593,1,-109.46795465403186,-335.09403393338613,0.01,0.01 -1581.6326530612246,1,-39.83805845141467,-276.40672308912804,0.01,0.01 -1632.6530612244899,1,24.75403165429632,-216.74299816205493,0.01,0.01 -1683.6734693877552,1,83.6775818280479,-157.9747940702846,0.01,0.01 -1734.6938775510205,1,136.8346222980739,-101.23656814885118,0.01,0.01 -1785.7142857142858,1,184.42422735042004,-47.16014775087199,0.01,0.01 -1836.734693877551,1,226.791523487859,3.949281188067366,0.01,0.01 -1887.7551020408164,1,264.33688157300367,51.99025505435129,0.01,0.01 -1938.7755102040817,1,297.4642986062226,96.98271751588322,0.01,0.01 -1989.795918367347,1,326.5539316633201,139.01437721880228,0.01,0.01 -2040.8163265306123,1,351.9489088146203,178.2064194774199,0.01,0.01 -2091.8367346938776,1,373.95024780709224,214.6919203425866,0.01,0.01 -2142.857142857143,1,392.8161509989211,248.60243215920502,0.01,0.01 -2193.877551020408,1,408.7634728691813,280.0597315825623,0.01,0.01 -2244.8979591836737,1,421.9700856579889,309.1707572133898,0.01,0.01 -2295.918367346939,1,432.57742176960176,336.0244485027215,0.01,0.01 -2346.9387755102043,1,440.69279444539876,360.68964087424973,0.01,0.01 -2397.9591836734694,1,446.39128232501747,383.2134537762328,0.01,0.01 -2448.979591836735,1,449.7170656626525,403.6197856857568,0.01,0.01 -2500.0,1,450.6841543184894,421.90763264506023,0.01,0.01 +0.0,1,-44.97595264191645,-40.40063509461096,0.01,0.01 +0.02040816326530612,1,-44.30820072074837,-45.71952453300766,0.01,0.01 +0.04081632653061224,1,-41.04700178820174,-48.201058902486146,0.01,0.01 +0.061224489795918366,1,-35.74054958085938,-47.965832008749146,0.01,0.01 +0.08163265306122448,1,-29.106223608047067,-45.48435880817074,0.01,0.01 +0.1020408163265306,1,-21.8253835127802,-41.370259974357374,0.01,0.01 +0.12244897959183673,1,-14.428317330209142,-36.211598246977516,0.01,0.01 +0.14285714285714285,1,-7.2699135408338,-30.48491432830629,0.01,0.01 +0.16326530612244897,1,-0.5565700744014579,-24.536472548444245,0.01,0.01 +0.18367346938775508,1,5.611711944213546,-18.598583371815057,0.01,0.01 +0.2040816326530612,1,11.203113265151215,-12.816822522196352,0.01,0.01 +0.22448979591836732,1,16.22563205150931,-7.27593756563277,0.01,0.01 +0.24489795918367346,1,20.708218518339898,-2.020291901969007,0.01,0.01 +0.26530612244897955,1,24.68910142063662,2.93143637707865,0.01,0.01 +0.2857142857142857,1,28.208969483025083,7.576234944513116,0.01,0.01 +0.3061224489795918,1,31.307217197702858,11.920150575624033,0.01,0.01 +0.32653061224489793,1,34.02003706808991,15.974029767980111,0.01,0.01 +0.3469387755102041,1,36.37958133646282,19.75081855840348,0.01,0.01 +0.36734693877551017,1,38.41371660681627,23.263874472541467,0.01,0.01 +0.3877551020408163,1,40.14608684277597,26.525923557689968,0.01,0.01 +0.4081632653061224,1,41.59631855303417,29.548420660047412,0.01,0.01 +0.42857142857142855,1,42.7802730988534,32.34115487431027,0.01,0.01 +0.44897959183673464,1,43.71029297512206,34.91199683952046,0.01,0.01 +0.4693877551020408,1,44.3954130505847,37.266719491329106,0.01,0.01 +0.4897959183673469,1,44.8415215664951,39.408846678587985,0.01,0.01 +0.5102040816326531,1,45.05146258834253,41.339496315275746,0.01,0.01 +0.5306122448979591,1,45.02507581579308,43.057194255088284,0.01,0.01 +0.5510204081632653,1,44.75917087078464,44.55763784084132,0.01,0.01 +0.5714285714285714,1,44.247434007528476,45.83339021912297,0.01,0.01 +0.5918367346938775,1,43.480265674963206,46.87348554962607,0.01,0.01 +0.6122448979591836,1,42.444548173623666,47.66292429386228,0.01,0.01 +0.6326530612244897,1,41.12334547730126,48.182034574487595,0.01,0.01 +0.6530612244897959,1,39.495542722455234,48.40567440725761,0.01,0.01 +0.673469387755102,1,37.53544417786146,48.302249377305294,0.01,0.01 +0.6938775510204082,1,35.212369445793755,47.832525882942925,0.01,0.01 +0.7142857142857142,1,32.49032544307804,46.948236929022904,0.01,0.01 +0.7346938775510203,1,29.32789804090339,45.5905171549874,0.01,0.01 +0.7551020408163265,1,25.678621015774937,43.68828706664445,0.01,0.01 +0.7755102040816326,1,21.49227022844291,41.156869864693746,0.01,0.01 +0.7959183673469387,1,16.717838261070998,37.89742910959155,0.01,0.01 +0.8163265306122448,1,11.309414286826271,33.79835583006353,0.01,0.01 +0.836734693877551,1,5.236842173178732,28.74062952787787,0.01,0.01 +0.8571428571428571,1,-1.4962572079392693,22.610509345316203,0.01,0.01 +0.8775510204081632,1,-8.824224576507511,15.324492907047599,0.01,0.01 +0.8979591836734693,1,-16.579434040144825,6.872290977814089,0.01,0.01 +0.9183673469387754,1,-24.445019923631108,-2.619219253000618,0.01,0.01 +0.9387755102040816,1,-31.917621523654393,-12.809602893650707,0.01,0.01 +0.9591836734693877,1,-38.3162586778345,-23.088729032293035,0.01,0.01 +0.9795918367346939,1,-42.88164071438466,-32.600661274349676,0.01,0.01 +1.0,1,-44.97595264191646,-40.40063509461095,0.01,0.01 From 7048e56f358a68cfee8f5cf5c1d4abd75362f63d Mon Sep 17 00:00:00 2001 From: Farrah Date: Fri, 6 Jun 2025 10:40:44 -0700 Subject: [PATCH 44/44] HDF5 for shorted pos data set --- HDF5_shortpos.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 HDF5_shortpos.py diff --git a/HDF5_shortpos.py b/HDF5_shortpos.py new file mode 100644 index 00000000..b78dca41 --- /dev/null +++ b/HDF5_shortpos.py @@ -0,0 +1,20 @@ +from orbitize import results +import matplotlib.pyplot as plt +hdf5_filename = "/home/fmolina/orbitize/tests/my_posterior_shortpos.hdf5" + + +loaded_results = results.Results() # Create blank results object for loading +loaded_results.load_results(hdf5_filename) + +print(loaded_results.labels) + + +orbit_plot_fig = loaded_results.plot_orbits( + object_to_plot=1, + num_orbits_to_plot=100, + start_mjd=0, + sep_pa_end_year=1859 +) +# save fig +orbit_plot_fig.savefig("/home/fmolina/orbitize/reloaded_shortpos_plot.png") +plt.show() \ No newline at end of file