diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index fde1386b..2bf46701 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.6","generation_timestamp":"2024-11-19T07:42:51","documenter_version":"1.8.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.6","generation_timestamp":"2024-11-19T19:20:07","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/HD82134-plot-grid.png b/dev/HD82134-plot-grid.png index aba91858..7e270a17 100644 Binary files a/dev/HD82134-plot-grid.png and b/dev/HD82134-plot-grid.png differ diff --git a/dev/HD91312_pma_rv_astrom-plot-grid.png b/dev/HD91312_pma_rv_astrom-plot-grid.png index e4d4bbaf..e76d33c3 100644 Binary files a/dev/HD91312_pma_rv_astrom-plot-grid.png and b/dev/HD91312_pma_rv_astrom-plot-grid.png differ diff --git a/dev/HR8799_res_co-plot-grid.png b/dev/HR8799_res_co-plot-grid.png index 8a846134..57039a7b 100644 Binary files a/dev/HR8799_res_co-plot-grid.png and b/dev/HR8799_res_co-plot-grid.png differ diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47909.0937mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47909.0937mjd.txt index 1c44b3f7..c6d94d66 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47909.0937mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47909.0937mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:14 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:32 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47927.319675mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47927.319675mjd.txt index 52f8aff8..ac75e873 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47927.319675mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47927.319675mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:15 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:33 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47955.736125mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47955.736125mjd.txt index 64b1265d..e00c7b33 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47955.736125mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47955.736125mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:15 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:34 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47994.817875mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47994.817875mjd.txt index b5816fae..7f3487fe 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47994.817875mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47994.817875mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:16 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:34 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47995.2927mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47995.2927mjd.txt index f81f168c..9d60f013 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47995.2927mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47995.2927mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:16 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:34 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47995.731mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47995.731mjd.txt index 6b49b187..de8093b4 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47995.731mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47995.731mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:16 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:34 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47996.1693mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47996.1693mjd.txt index ccd29c4a..c6638554 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47996.1693mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47996.1693mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:16 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:34 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47996.6076mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47996.6076mjd.txt index 3cb3b9d8..507b7622 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47996.6076mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47996.6076mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:16 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:35 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-47997.0459mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-47997.0459mjd.txt index 8599b657..96396494 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-47997.0459mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-47997.0459mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:17 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:35 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48092.047425mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48092.047425mjd.txt index e8e6044d..60bf5f8d 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48092.047425mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48092.047425mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:17 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:35 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48092.485725mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48092.485725mjd.txt index e70cd791..c69069d0 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48092.485725mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48092.485725mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:17 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:35 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48106.69395mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48106.69395mjd.txt index c5973407..d5575693 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48106.69395mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48106.69395mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:17 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:35 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48139.56645mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48139.56645mjd.txt index 04f01790..cbd5b4bf 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48139.56645mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48139.56645mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:17 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:36 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48177.735075mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48177.735075mjd.txt index 8c199aad..632bff1e 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48177.735075mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48177.735075mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:17 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:36 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48178.173375mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48178.173375mjd.txt index 23aafacc..2a7d0921 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48178.173375mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48178.173375mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:18 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:36 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48178.611675mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48178.611675mjd.txt index c2ce52b3..f7eca203 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48178.611675mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48178.611675mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:18 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:36 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48182.629425mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48182.629425mjd.txt index d5ab17d3..f7fd07a6 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48182.629425mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48182.629425mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:18 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:36 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48183.067725mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48183.067725mjd.txt index 33f76778..a7d3af01 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48183.067725mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48183.067725mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:18 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:37 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48183.506025mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48183.506025mjd.txt index 2b217e02..99989c00 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48183.506025mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48183.506025mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:18 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:37 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48301.993125mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48301.993125mjd.txt index bd5092a4..f58a91a7 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48301.993125mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48301.993125mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:19 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:37 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48331.724475mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48331.724475mjd.txt index fa3b7802..564a237d 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48331.724475mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48331.724475mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:19 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:37 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48348.599025mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48348.599025mjd.txt index 840bb910..3dfb6dc8 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48348.599025mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48348.599025mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:19 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:37 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48485.275575mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48485.275575mjd.txt index 9b70b6e6..8ccfb99b 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48485.275575mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48485.275575mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:19 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:38 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48512.450175mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48512.450175mjd.txt index 5a714d25..83b69163 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48512.450175mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48512.450175mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:19 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:38 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48648.57885mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48648.57885mjd.txt index 71f9ee7b..ab66093e 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48648.57885mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48648.57885mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:20 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:38 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48671.1513mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48671.1513mjd.txt index 6bc4c567..7c33d2cd 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48671.1513mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48671.1513mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:20 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:38 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48695.55mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48695.55mjd.txt index ecfe6955..b2a2eba5 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48695.55mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48695.55mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:20 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:38 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/_geocentre_pos/HORIZONS-Geocenter-48831.7152mjd.txt b/dev/_geocentre_pos/HORIZONS-Geocenter-48831.7152mjd.txt index 25c65a8d..4289e03c 100644 --- a/dev/_geocentre_pos/HORIZONS-Geocenter-48831.7152mjd.txt +++ b/dev/_geocentre_pos/HORIZONS-Geocenter-48831.7152mjd.txt @@ -30,7 +30,7 @@ ******************************************************************************* -Ephemeris / API_USER Mon Nov 18 21:56:20 2024 Pasadena, USA / Horizons +Ephemeris / API_USER Tue Nov 19 09:32:39 2024 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE441} Center body name: Solar System Barycenter (0) {source: DE441} diff --git a/dev/api/index.html b/dev/api/index.html index e3b9d715..341d0de2 100644 --- a/dev/api/index.html +++ b/dev/api/index.html @@ -3,15 +3,15 @@ [prior_1] ~ [UnivariateDistribution] [prior_2] ~ [UnivariateDistribution] calculation_3 = [planet_name].[prior_1] + [planet_name].[prior_2] + [system.variable_x] -end [likelihood_objects...]

Generate a Planet model named planet_name. A variable will be created with the name [planet_name] in the current scope. Orbit_Type specifies the orbit parameterization from PlanetOrbits.jl. You must provide all input variables needed for the selected orbit type (see PlanetOrbits.jl documentation). Following that is a block of variable assignments. Variables with a ~ will be free variables with a prior distribution given by the right-hand-side (a UnivariateDistribution from Distributions.jl or a KDEDist). Calculated quantities are also allowed. These may reference other variables using the planet name followed by a dot and the name of the variable. Variables from other planets in a single system are not accessible. You can access other variables in the current local scope, but these bindings are only guaranteed to be resolved a single time. Note that using non-constant global variables in calculated expressions can lead to poor performance. Finally, the planet model can be conditioned on data by supplying zero or more likelihood objects.

source
Octofitter.@systemMacro
@system [system_name] begin
+end [likelihood_objects...]

Generate a Planet model named planet_name. A variable will be created with the name [planet_name] in the current scope. Orbit_Type specifies the orbit parameterization from PlanetOrbits.jl. You must provide all input variables needed for the selected orbit type (see PlanetOrbits.jl documentation). Following that is a block of variable assignments. Variables with a ~ will be free variables with a prior distribution given by the right-hand-side (a UnivariateDistribution from Distributions.jl or a KDEDist). Calculated quantities are also allowed. These may reference other variables using the planet name followed by a dot and the name of the variable. Variables from other planets in a single system are not accessible. You can access other variables in the current local scope, but these bindings are only guaranteed to be resolved a single time. Note that using non-constant global variables in calculated expressions can lead to poor performance. Finally, the planet model can be conditioned on data by supplying zero or more likelihood objects.

source
Octofitter.@systemMacro
@system [system_name] begin
     [prior_1] ~ [UnivariateDistribution]
     [prior_2] ~ [UnivariateDistribution]
     calculation_3 = system.[prior_1] + system.[prior_2]
-end [likelihood_objects...] [planet_models...]

Generate a System model named system_name. A variable will be created with the name [system_name] in the current scope. Following that is a block of variable assignments. Variables with a ~ will be free variables with a prior distribution given by the right-hand-side (a UnivariateDistribution from Distributions.jl or a KDEDist). Calculated quantities are also allowed. These may reference other variables using the planet name followed by a dot and the name of the variable. Variables from other planets in a single system are not accessible. You can access other variables in the current local scope, but these bindings are only guaranteed to be resolved a single time. Note that using non-constant global variables in calculated expressions can lead to poor performance. After the end of the variable block, the system model can be conditioned on data by supplying zero or more likelihood objects. Finally, zero or more planet models can be attached to the system, potentially conditioned on likelihood objects of their own.

source
Octofitter.PlanetRelAstromLikelihoodType
PlanetRelAstromLikelihood(
+end [likelihood_objects...] [planet_models...]

Generate a System model named system_name. A variable will be created with the name [system_name] in the current scope. Following that is a block of variable assignments. Variables with a ~ will be free variables with a prior distribution given by the right-hand-side (a UnivariateDistribution from Distributions.jl or a KDEDist). Calculated quantities are also allowed. These may reference other variables using the planet name followed by a dot and the name of the variable. Variables from other planets in a single system are not accessible. You can access other variables in the current local scope, but these bindings are only guaranteed to be resolved a single time. Note that using non-constant global variables in calculated expressions can lead to poor performance. After the end of the variable block, the system model can be conditioned on data by supplying zero or more likelihood objects. Finally, zero or more planet models can be attached to the system, potentially conditioned on likelihood objects of their own.

source
Octofitter.PlanetRelAstromLikelihoodType
PlanetRelAstromLikelihood(
     (epoch = 5000, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
     (epoch = 5050, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
     (epoch = 5100, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
-)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch is a required column, in addition to either :ra, :dec, :σ_ra, :σ_dec or :pa, :sep, :σ_pa, :σ_sep. All units are in milliarcseconds or radians as appropriate.

In addition to the example above, any Tables.jl compatible source can be provided.

source
OctofitterRadialVelocity.StarAbsoluteRVLikelihoodType
StarAbsoluteRVLikelihood(
+)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch is a required column, in addition to either :ra, :dec, :σ_ra, :σ_dec or :pa, :sep, :σ_pa, :σ_sep. All units are in milliarcseconds or radians as appropriate.

In addition to the example above, any Tables.jl compatible source can be provided.

source
OctofitterRadialVelocity.StarAbsoluteRVLikelihoodType
StarAbsoluteRVLikelihood(
     (;epoch=5000.0,  rv=−6.54, σ_rv=1.30),
     (;epoch=5050.1,  rv=−3.33, σ_rv=1.09),
     (;epoch=5100.2,  rv=7.90,  σ_rv=.11),
@@ -23,18 +23,18 @@
     # Optional:
     trend_function=(θ_system, epoch)->0.,
     gaussian_process=nothing
-)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch (mjd), :rv (m/s), and :σ_rv (m/s) are all required.

The offset and jitter parameters specify which variables should be read from the model for the RV zero-point and jitter of this instrument.

In addition to the example above, any Tables.jl compatible source can be provided.

source
OctofitterRadialVelocity.PlanetRelativeRVLikelihoodType
PlanetRelativeRVLikelihood(
+)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch (mjd), :rv (m/s), and :σ_rv (m/s) are all required.

The offset and jitter parameters specify which variables should be read from the model for the RV zero-point and jitter of this instrument.

In addition to the example above, any Tables.jl compatible source can be provided.

source
OctofitterRadialVelocity.PlanetRelativeRVLikelihoodType
PlanetRelativeRVLikelihood(
     (;epoch=5000.0,  rv=−6.54, σ_rv=1.30),
     (;epoch=5050.1,  rv=−3.33, σ_rv=1.09),
     (;epoch=5100.2,  rv=7.90,  σ_rv=.11),
 
     jitter=:jitter_1,
     instrument_name="inst name",
-)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch (mjd), :rv (m/s), and :σ_rv (m/s) are all required.

In addition to the example above, any Tables.jl compatible source can be provided.

The jitter parameter specify which variables should be read from the model for the jitter of this instrument.

source
Octofitter.PhotometryLikelihoodType
PhotometryLikelihood(
+)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch (mjd), :rv (m/s), and :σ_rv (m/s) are all required.

In addition to the example above, any Tables.jl compatible source can be provided.

The jitter parameter specify which variables should be read from the model for the jitter of this instrument.

source
Octofitter.PhotometryLikelihoodType
PhotometryLikelihood(
     (band = :Z, phot=15.0, σ_phot=3.),
     (band = :J, phot=13.5, σ_phot=0.5),
     (band = :L, phot=11.0, σ_phot=1.0)
-)

A likelihood for comparing measured photometry points in one or more filter bands to data (provided here). Requires the :band, :phot', and:σ_phot` columns. Can be provided with any Tables.jl compatible data source.

source
Octofitter.HGCALikelihoodType
HGCALikelihood(;gaia_id=1234)

Load proper motion anomaly data from the HIPPARCOS-GAIA Catalog of Accelerations (Brandt et al) for a star with catalog id gaia_id. The resulting velocities are in mas/yr and have the long term trend between HIPPARCOS and GAIA already subtracted out. e.g. we would expect 0 pma if there is no companion.

source
Octofitter.ObsPriorAstromONeil2019Type
ObsPriorAstromONeil2019(astrometry_likelihood, period_prior)

Given a an astrometry likelihood (PlanetRelAstromLikelihood), apply the "observable based priors" of K. O'Neil 2019 "Improving Orbit Estimates for Incomplete Orbits with a New Approach to Priors: with Applications from Black Holes to Planets".

This prior correction is only correct if you supply Uniform priors on all Campbell orbital parameters and a Uniform prior on Period (not semi-major axis). This period prior has a significant impact in the fit and recommendations for its range were not published in the original paper.

Examples

astrom_like = PlanetRelAstromLikelihood(astrom_table)
+)

A likelihood for comparing measured photometry points in one or more filter bands to data (provided here). Requires the :band, :phot', and:σ_phot` columns. Can be provided with any Tables.jl compatible data source.

source
Octofitter.HGCALikelihoodType
HGCALikelihood(;gaia_id=1234)

Load proper motion anomaly data from the HIPPARCOS-GAIA Catalog of Accelerations (Brandt et al) for a star with catalog id gaia_id. The resulting velocities are in mas/yr and have the long term trend between HIPPARCOS and GAIA already subtracted out. e.g. we would expect 0 pma if there is no companion.

source
Octofitter.ObsPriorAstromONeil2019Type
ObsPriorAstromONeil2019(astrometry_likelihood, period_prior)

Given a an astrometry likelihood (PlanetRelAstromLikelihood), apply the "observable based priors" of K. O'Neil 2019 "Improving Orbit Estimates for Incomplete Orbits with a New Approach to Priors: with Applications from Black Holes to Planets".

This prior correction is only correct if you supply Uniform priors on all Campbell orbital parameters and a Uniform prior on Period (not semi-major axis). This period prior has a significant impact in the fit and recommendations for its range were not published in the original paper.

Examples

astrom_like = PlanetRelAstromLikelihood(astrom_table)
 
 # Apply observable based priors ontop of our uniform Campbell priors:
 obs_prior = ObsPriorAstromONeil2019(astrom_like)
@@ -59,7 +59,7 @@
     ω ~ UniformCircular()
     Ω ~ UniformCircular()
     τ ~ UniformCircular(1.0)
-end astrom_like obs_prior
source
Octofitter.SineType
Sine()

A custom univariate distribution. The pdf is a sine function defined between 0 and π. This is a common prior distribution used when fitting orbits to astrometry.

The full Distributions.jl interface is not yet defined for this distribution, but the following methods work: pdf, logpdf, minimum, maximum, insupport, mean, var, cdf, quantile

source
PlanetOrbits.mjdFunction
mjd("2020-01-01")

Get the modfied julian day of a date, or in general a UTC timestamp.

source
mjd(Date("2020-01-01"))

Get the modfied julian day of a Date or DateTime object.

source
mjd()

Get the current modified julian day of right now.

source
Octofitter.gaia_plxFunction
gaia_plx(gaia_id=12123)

Get a distribution (truncated Normal) of parallax distance in mas of a source with GAIA catalog id gaia_id.

source
Octofitter.advancedhmcFunction

The method signature of Octofitter.hmc is as follows:

advancedhmc(
+end astrom_like obs_prior
source
Octofitter.SineType
Sine()

A custom univariate distribution. The pdf is a sine function defined between 0 and π. This is a common prior distribution used when fitting orbits to astrometry.

The full Distributions.jl interface is not yet defined for this distribution, but the following methods work: pdf, logpdf, minimum, maximum, insupport, mean, var, cdf, quantile

source
PlanetOrbits.mjdFunction
mjd("2020-01-01")

Get the modfied julian day of a date, or in general a UTC timestamp.

source
mjd(Date("2020-01-01"))

Get the modfied julian day of a Date or DateTime object.

source
mjd()

Get the current modified julian day of right now.

source
Octofitter.gaia_plxFunction
gaia_plx(gaia_id=12123)

Get a distribution (truncated Normal) of parallax distance in mas of a source with GAIA catalog id gaia_id.

source
Octofitter.advancedhmcFunction

The method signature of Octofitter.hmc is as follows:

advancedhmc(
     [rng::Random.AbstractRNG],
     model::Octofitter.LogDensityModel
     target_accept::Number=0.8,
@@ -72,7 +72,7 @@
     initial_parameters=nothing, # deprecated
     step_size=nothing,
     verbosity=2,
-)

The only required arguments are system, adaptation, and iterations. The two positional arguments are system, the model you wish to sample; and targetaccept, the acceptance rate that should be targeted during windowed adaptation. During this time, the step size and mass matrix will be adapted (see AdvancedHMC.jl for more information). The number of steps taken during adaptation is controlled by adaptation. You can prevent these samples from being dropped by pasing includeadaptation=false. The total number of posterior samples produced are given by iterations. These include the adaptation steps that may be discarded. treedepth controls the maximum tree depth of the sampler. initialparameters is an optional way to pass a starting point for the chain. If you don't pass a default position, one will be selected by drawing initial_samples from the priors. The sample with the highest posterior value will be used as the starting point.

source
PlanetOrbits.VisualType
Visual{OrbitType}(..., plx=...)

This wraps another orbit to add the parallax distance field plx, thus allowing projected quantities to be calculated. It forwards everything else to the parent orbit.

For example, the KepOrbit type supports calculating x and y positions in AU. A Visual{KepOrbit} additionally supports calculating projected right ascension and declination offsets.

Note

The ThieleInnesOrbit type does not need to be wrapped in Visual as it the Thiele-Innes constants are already expressed in milliarcseconds and thus it always requires a plx value.

source
Octofitter.sonora_photometry_interpolatorFunction
sonora_photometry_interpolator(:Keck_L′, [metalicity="+0.0"])

Given a supported photometric band and [M/H] metalicity (default=solar), return a function of temperature (K) and mass (M_jup) that gives the absolute magnitude of the planet in that bandpass.

Supported bands: :MKOY, :MKOZ, :MKOJ, :MKOH, :MKOK, :MKOL′, :MKOM′, :TwoMASSJ, :TwoMASSH, :TwoMASSKs, :KeckKs, :KeckL′, :KeckMs, :SDSSg′, :SDSSr′, :SDSSi′, :SDSSz′, :IRAC36, :IRAC45, :IRAC57, :IRAC79, :WISEW1, :WISEW2, :WISEW3, :WISE_W4

Supported metalicities: "+0.0", "-0.5", "+0.5"

source
Octofitter.sonora_cooling_interpolatorFunction
itp = sonora_cooling_interpolator()

Create a function mapping (ageMyr, massMjup) -> temp_K using Sonora Bobcat cooling model grids.

source
Octofitter.plotchainsFunction
plotchains(
+)

The only required arguments are system, adaptation, and iterations. The two positional arguments are system, the model you wish to sample; and targetaccept, the acceptance rate that should be targeted during windowed adaptation. During this time, the step size and mass matrix will be adapted (see AdvancedHMC.jl for more information). The number of steps taken during adaptation is controlled by adaptation. You can prevent these samples from being dropped by pasing includeadaptation=false. The total number of posterior samples produced are given by iterations. These include the adaptation steps that may be discarded. treedepth controls the maximum tree depth of the sampler. initialparameters is an optional way to pass a starting point for the chain. If you don't pass a default position, one will be selected by drawing initial_samples from the priors. The sample with the highest posterior value will be used as the starting point.

source
PlanetOrbits.VisualType
Visual{OrbitType}(..., plx=...)

This wraps another orbit to add the parallax distance field plx, thus allowing projected quantities to be calculated. It forwards everything else to the parent orbit.

For example, the KepOrbit type supports calculating x and y positions in AU. A Visual{KepOrbit} additionally supports calculating projected right ascension and declination offsets.

Note

The ThieleInnesOrbit type does not need to be wrapped in Visual as it the Thiele-Innes constants are already expressed in milliarcseconds and thus it always requires a plx value.

source
Octofitter.sonora_photometry_interpolatorFunction
sonora_photometry_interpolator(:Keck_L′, [metalicity="+0.0"])

Given a supported photometric band and [M/H] metalicity (default=solar), return a function of temperature (K) and mass (M_jup) that gives the absolute magnitude of the planet in that bandpass.

Supported bands: :MKOY, :MKOZ, :MKOJ, :MKOH, :MKOK, :MKOL′, :MKOM′, :TwoMASSJ, :TwoMASSH, :TwoMASSKs, :KeckKs, :KeckL′, :KeckMs, :SDSSg′, :SDSSr′, :SDSSi′, :SDSSz′, :IRAC36, :IRAC45, :IRAC57, :IRAC79, :WISEW1, :WISEW2, :WISEW3, :WISE_W4

Supported metalicities: "+0.0", "-0.5", "+0.5"

source
Octofitter.sonora_cooling_interpolatorFunction
itp = sonora_cooling_interpolator()

Create a function mapping (ageMyr, massMjup) -> temp_K using Sonora Bobcat cooling model grids.

source
Octofitter.plotchainsFunction
plotchains(
     chain, planet_key;
     N=1500,
     ii = rand(1:size(chain,1)*size(chain,3), N),
@@ -83,4 +83,4 @@
     alpha=30/length(ii),
     attime=nothing,
     kwargs...,
-)

Draw samples from a posterior chain for a given planet given by name planet_key and visualize them in some way. Use kind to control what plot is made. A few options: :astrometry, :radvel, :trueanom, :meananom, :eccanom, :x, :y, :z, (:x, :y), :raoff, :decoff, :pmra, :pmdec, :accra, :accdec, :radvel, :posangle, :projectedseparation. See PlanetOrbits documentation for more details.

Inputs:

  • chain The chain to draw from
  • planet_key Planet name in the model (symbol)
  • N=1500 Number of samples to draw for the plot
  • kind=nothing Specify what kind of plot to make.
  • ii=... Specific row numbers to use, if you want to e.g. plot the same 100 samples in a few different plots
  • color="planetkeya" Column name to to map colors to. Semi-major axis by default but can be any column or an arbitrary array.
  • colorbartitle=color Name for colourbar
  • clims=nothing Tuple of colour limits (min and max)
  • cmap=:plasma Colormap
  • alpha=... Transparency of the lines
source
Octofitter.projectpositionsFunction
projectpositions(chains.planets[1], mjd("2020-02-02"))

Given the posterior for a particular planet in the model and a modified julian date(s), return ra and dec offsets in mas for each sampling in the posterior.

source
+)

Draw samples from a posterior chain for a given planet given by name planet_key and visualize them in some way. Use kind to control what plot is made. A few options: :astrometry, :radvel, :trueanom, :meananom, :eccanom, :x, :y, :z, (:x, :y), :raoff, :decoff, :pmra, :pmdec, :accra, :accdec, :radvel, :posangle, :projectedseparation. See PlanetOrbits documentation for more details.

Inputs:

source
Octofitter.projectpositionsFunction
projectpositions(chains.planets[1], mjd("2020-02-02"))

Given the posterior for a particular planet in the model and a modified julian date(s), return ra and dec offsets in mas for each sampling in the posterior.

source
diff --git a/dev/astrom-pma-rv/08833924.png b/dev/astrom-pma-rv/08833924.png new file mode 100644 index 00000000..e76d33c3 Binary files /dev/null and b/dev/astrom-pma-rv/08833924.png differ diff --git a/dev/astrom-pma-rv/6406e6d8.png b/dev/astrom-pma-rv/6406e6d8.png deleted file mode 100644 index e4d4bbaf..00000000 Binary files a/dev/astrom-pma-rv/6406e6d8.png and /dev/null differ diff --git a/dev/astrom-pma-rv/index.html b/dev/astrom-pma-rv/index.html index a775dadb..8e17a51a 100644 --- a/dev/astrom-pma-rv/index.html +++ b/dev/astrom-pma-rv/index.html @@ -50,7 +50,7 @@ model_pma = Octofitter.LogDensityModel(HD91312_pma)
LogDensityModel for System HD91312_pma of dimension 14 and 7 epochs with fields .ℓπcallback and .∇ℓπcallback
 

After the priors, we add the proper motion anomaly measurements from the HGCA. If this is your first time running this code, you will be prompted to automatically download and cache the catalog which may take around 30 seconds.

Sampling from the posterior (PMA only)

Because proper motion anomaly data is quite sparse, it can often produce multi-modal posteriors. If your orbit already has several relative astrometry or RV data points, this is less of an issue. But in many cases it is recommended to use the Pigeons.jl sampler instead of Octofitter's default. This sampler is less efficient for unimodal distributions, but is more robust at exploring posteriors with distinct, widely separated peaks.

To install and use Pigeons.jl with Octofitter, type using Pigeons at in the terminal and accept the prompt to install the package. You may have to restart Julia.

Note

octofit_pigeons scales very well across multiple cores. Start julia with julia --threads=auto to make sure you have multiple threads available for sampling.

We now sample from our model using Pigeons:

using Pigeons
 chain_pma, pt = octofit_pigeons(model_pma, n_rounds=13, explorer=SliceSampler())
-display(chain_pma)
┌ Warning: Module Octofitter with build ID fafbfcfd-c563-0173-0000-003d53fb7590 is missing from the cache.
+display(chain_pma)
┌ Warning: Module Octofitter with build ID fafbfcfd-1710-b70a-0000-003b9bab10ca is missing from the cache.
 │ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
 └ @ Base loading.jl:1948
 ┌ Warning: Module ForwardDiff with build ID fafbfcfd-80d5-a3a4-0000-0054f268105a is missing from the cache.
@@ -64,19 +64,19 @@
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
 ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
-        2          0       4.53        2.1       1.04   2.42e+07  -2.55e+03          0      0.786          1          1
-        4          0       4.51       6.93      0.141    1.5e+05      -20.9          0      0.631          1          1
+        2          0       4.53        2.1      0.752   2.42e+07  -2.55e+03          0      0.786          1          1
+        4          0       4.51       6.93      0.118    1.5e+05      -20.9          0      0.631          1          1
         8          0       6.07       7.48      0.225   2.53e+05      -15.5   2.52e-75      0.563          1          1
-       16          0       6.51       8.29      0.434   3.76e+05      -17.6   1.53e-45      0.523          1          1
-       32          0       7.52       9.24      0.853   6.21e+05      -20.1   9.87e-08      0.459          1          1
-       64          0       7.68       6.44       1.89   6.15e+07      -24.6     0.0685      0.544          1          1
-      128          2       8.63       6.65        3.3   1.06e+08      -24.6      0.169      0.507          1          1
-      256          4       8.83       6.45       6.59   2.13e+08      -22.6       0.35      0.507          1          1
-      512         20       8.67       6.49       13.2   4.23e+08      -23.2      0.327      0.511          1          1
- 1.02e+03         51       8.78       6.23       26.4   8.55e+08      -22.5       0.37      0.516          1          1
- 2.05e+03         99       8.66       6.01       53.3    1.7e+09        -23      0.398      0.527          1          1
-  4.1e+03        197       8.74       6.18        106   3.42e+09      -22.5      0.394      0.519          1          1
- 8.19e+03        415       8.71       6.19        213   6.82e+09      -22.5      0.409      0.519          1          1
+       16          0       6.51       8.29      0.436   3.76e+05      -17.6   1.53e-45      0.523          1          1
+       32          0       7.52       9.24      0.858   6.21e+05      -20.1   9.87e-08      0.459          1          1
+       64          0       7.68       6.44       1.93   6.15e+07      -24.6     0.0685      0.544          1          1
+      128          2       8.63       6.65       3.35   1.06e+08      -24.6      0.169      0.507          1          1
+      256          4       8.83       6.45        6.7   2.13e+08      -22.6       0.35      0.507          1          1
+      512         20       8.67       6.49       13.3   4.23e+08      -23.2      0.327      0.511          1          1
+ 1.02e+03         51       8.78       6.23       26.7   8.55e+08      -22.5       0.37      0.516          1          1
+ 2.05e+03         99       8.66       6.01       53.7    1.7e+09        -23      0.398      0.527          1          1
+  4.1e+03        197       8.74       6.18        107   3.42e+09      -22.5      0.394      0.519          1          1
+ 8.19e+03        415       8.71       6.19        215   6.82e+09      -22.5      0.409      0.519          1          1
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Note that octofit_pigeons took somewhat longer to run than octofit typically does; however, as we will see, it sampled successfully from severally completely disconnected modes in the posterior. That makes it a good fit for sampling from proper motion anomaly and relative astrometry with limited orbital coverage.

Pair Plot

If we wish to examine the covariance between parameters in more detail, we can construct a pair-plot (aka. corner plot).

# Create a corner plot / pair plot.
 # We can access any property from the chain specified in Variables
 using CairoMakie
@@ -134,13 +134,13 @@
 chain_pma_astrom, pt = octofit_pigeons(model_pma_astrom, n_rounds=12, explorer=SliceSampler())
[ Info: Preparing model
 ┌ Info: Determined number of free variables
 └   D = 14
-ℓπcallback(θ, args...): 0.000018 seconds (9 allocations: 1.688 KiB)
+ℓπcallback(θ, args...): 0.000032 seconds (9 allocations: 1.688 KiB)
 ┌ Info: Tuning autodiff
 │   chunk_size = 14
 └   t = 1.2774e-5
 ┌ Info: Selected auto-diff chunk size
 └   ideal_chunk_size = 14
-∇ℓπcallback(θ): 0.000050 seconds (1 allocation: 176 bytes)
+∇ℓπcallback(θ): 0.000051 seconds (1 allocation: 176 bytes)
 [ Info: Sampler running with multiple threads     : true
 [ Info: Likelihood evaluated with multiple threads: false
 [ Info: Determining initial positions and metric using pathfinder
@@ -149,18 +149,18 @@
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
 ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
-        2          0       3.61        4.2      0.321   4.02e+06  -6.54e+03          0      0.748          1          1
-        4          0       5.38       4.82      0.235   1.52e+05       -462          0      0.671          1          1
-        8          0        6.8       6.46      0.435   2.44e+05       -391   1.8e-292      0.572          1          1
-       16          0       6.89       7.44      0.843   4.01e+05      -99.7   1.58e-96      0.538          1          1
+        2          0       3.61        4.2      0.327   4.02e+06  -6.54e+03          0      0.748          1          1
+        4          0       5.38       4.82      0.236   1.52e+05       -462          0      0.671          1          1
+        8          0        6.8       6.46      0.443   2.44e+05       -391   1.8e-292      0.572          1          1
+       16          0       6.89       7.44      0.863   4.01e+05      -99.7   1.58e-96      0.538          1          1
        32          0       8.13        8.1       1.63   5.64e+05      -69.7   7.27e-19      0.477          1          1
-       64          1       8.66       4.96       3.41   6.57e+07      -76.8     0.0155      0.561          1          1
-      128          4       9.13       4.72       6.53   1.19e+08      -75.7    0.00247      0.553          1          1
+       64          1       8.66       4.96       3.44   6.57e+07      -76.8     0.0155      0.561          1          1
+      128          4       9.13       4.72       6.54   1.19e+08      -75.7    0.00247      0.553          1          1
       256          8       9.81       5.05       13.1   2.38e+08      -75.5     0.0166       0.52          1          1
-      512         15       9.96       6.68       25.7   4.66e+08      -75.1     0.0478      0.463          1          1
- 1.02e+03         22         10       6.49       51.3   9.22e+08      -74.1      0.127      0.468          1          1
+      512         15       9.96       6.68       25.8   4.66e+08      -75.1     0.0478      0.463          1          1
+ 1.02e+03         22         10       6.49       51.5   9.22e+08      -74.1      0.127      0.468          1          1
  2.05e+03         52       10.2       6.95        103   1.85e+09      -73.8       0.24      0.448          1          1
-  4.1e+03        107       10.1       6.96        203   3.64e+09      -73.6      0.292      0.449          1          1
+  4.1e+03        107       10.1       6.96        204   3.64e+09      -73.6      0.292      0.449          1          1
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Model: PMA & Relative Astrometry & RVs

We now add in three additional epochs of stellar RVs.

using OctofitterRadialVelocity
 
 rvlike = StarAbsoluteRVLikelihood(
@@ -207,13 +207,13 @@
 display(chain_pma_rv_astrom)
[ Info: Preparing model
 ┌ Info: Determined number of free variables
 └   D = 16
-ℓπcallback(θ, args...): 0.000021 seconds (9 allocations: 1.828 KiB)
+ℓπcallback(θ, args...): 0.000033 seconds (9 allocations: 1.828 KiB)
 ┌ Info: Tuning autodiff
 │   chunk_size = 16
-└   t = 1.6241e-5
+└   t = 1.638e-5
 ┌ Info: Selected auto-diff chunk size
 └   ideal_chunk_size = 16
-∇ℓπcallback(θ): 0.000057 seconds (1 allocation: 192 bytes)
+∇ℓπcallback(θ): 0.000082 seconds (1 allocation: 192 bytes)
 [ Info: Sampler running with multiple threads     : true
 [ Info: Likelihood evaluated with multiple threads: false
 [ Info: Determining initial positions and metric using pathfinder
@@ -222,18 +222,18 @@
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
 ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
-        2          0       3.56       3.63      0.408   6.98e+06  -2.85e+03          0      0.768          1          1
-        4          0       7.04       5.45      0.351   1.55e+05  -1.22e+03          0      0.597          1          1
-        8          0        5.4       7.23      0.675   2.49e+05       -428          0      0.592          1          1
-       16          0       8.06       6.86       1.29   3.97e+05       -122    4.1e-23      0.519          1          1
-       32          0       8.59       8.53       2.49   5.68e+05       -172   9.86e-66      0.448          1          1
-       64          0       10.2       4.94       5.21   7.81e+07        176   0.000942      0.511          1          1
-      128          0       9.62       5.62       10.1   1.44e+08       -117    0.00641      0.508          1          1
-      256          1       10.1       7.11       19.8   2.79e+08       -115    0.00863      0.445          1          1
-      512          4       10.2       7.78       39.2    5.5e+08       -113     0.0186       0.42          1          1
- 1.02e+03          8       10.6       8.08       77.6    1.1e+09       -108      0.044      0.399          1          1
- 2.05e+03         42       10.7       7.06        153   2.19e+09      -94.3     0.0949      0.428          1          1
-  4.1e+03         70       10.6       7.68        300    4.3e+09       -100      0.174      0.411          1          1
+        2          0       3.56       3.63      0.445   6.98e+06  -2.85e+03          0      0.768          1          1
+        4          0       7.04       5.45      0.359   1.55e+05  -1.22e+03          0      0.597          1          1
+        8          0        5.4       7.23      0.677   2.49e+05       -428          0      0.592          1          1
+       16          0       8.06       6.86        1.3   3.97e+05       -122    4.1e-23      0.519          1          1
+       32          0       8.59       8.53        2.5   5.68e+05       -172   9.86e-66      0.448          1          1
+       64          0       10.2       4.94       5.25   7.81e+07        176   0.000942      0.511          1          1
+      128          0       9.62       5.62       10.2   1.44e+08       -117    0.00641      0.508          1          1
+      256          1       10.1       7.11         20   2.79e+08       -115    0.00863      0.445          1          1
+      512          4       10.2       7.78       39.5    5.5e+08       -113     0.0186       0.42          1          1
+ 1.02e+03          8       10.6       8.08       78.2    1.1e+09       -108      0.044      0.399          1          1
+ 2.05e+03         42       10.7       7.06        154   2.19e+09      -94.3     0.0949      0.428          1          1
+  4.1e+03         70       10.6       7.68        302    4.3e+09       -100      0.174      0.411          1          1
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

The mass vs. semi-major axis posterior is now much more constrained:

using CairoMakie, PairPlots
 pairplot(
     (; a=chain_pma_rv_astrom["b_a"][:], mass=chain_pma_rv_astrom["b_mass"][:]) =>
@@ -243,7 +243,7 @@
             PairPlots.MarginConfidenceLimits()
         ),
     labels=Dict(:mass=>"mass [Mⱼᵤₚ]", :a=>"sma. [au]"),
-)
Example block output

It is now useful to display the orbits projected onto the plane of the sky using octoplot. This function produces a nine-panel figure showing posterior predictive distributions for velocity (in three dimensions), projected positions vs. time in the plane of the sky, and various other two and three-dimensional views.

octoplot(model_pma_rv_astrom, chain_pma_rv_astrom, show_mass=true)
Example block output

Model: Relative Astrometry & RVs (no PMA)

There is a final model we should consider: one using the RV and astrometry data, but not the proper motion anomaly:

using OctofitterRadialVelocity
+)
Example block output

It is now useful to display the orbits projected onto the plane of the sky using octoplot. This function produces a nine-panel figure showing posterior predictive distributions for velocity (in three dimensions), projected positions vs. time in the plane of the sky, and various other two and three-dimensional views.

octoplot(model_pma_rv_astrom, chain_pma_rv_astrom, show_mass=true)
Example block output

Model: Relative Astrometry & RVs (no PMA)

There is a final model we should consider: one using the RV and astrometry data, but not the proper motion anomaly:

using OctofitterRadialVelocity
 
 @planet b Visual{KepOrbit} begin
     a ~ LogUniform(0.1,400)
@@ -275,13 +275,13 @@
 chain_rv_astrom, pt = octofit_pigeons(model_rv_astrom, n_rounds=12)
[ Info: Preparing model
 ┌ Info: Determined number of free variables
 └   D = 14
-ℓπcallback(θ, args...): 0.000023 seconds (9 allocations: 608 bytes)
+ℓπcallback(θ, args...): 0.000034 seconds (9 allocations: 608 bytes)
 ┌ Info: Tuning autodiff
 │   chunk_size = 14
-└   t = 1.1461e-5
+└   t = 1.1342e-5
 ┌ Info: Selected auto-diff chunk size
 └   ideal_chunk_size = 14
-∇ℓπcallback(θ): 0.000041 seconds (1 allocation: 176 bytes)
+∇ℓπcallback(θ): 0.000053 seconds (1 allocation: 176 bytes)
 [ Info: Sampler running with multiple threads     : true
 [ Info: Likelihood evaluated with multiple threads: false
 [ Info: Determining initial positions and metric using pathfinder
@@ -290,18 +290,18 @@
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
 ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
-        2          0       4.56        3.8      0.308   4.03e+06       -216          0       0.73          1          1
-        4          0       3.75       3.29      0.215   1.63e+05       -108    1.8e-64      0.773          1          1
-        8          0       5.61       6.73       0.41    2.3e+05      -83.1   0.000173      0.602          1          1
-       16          0       6.78       7.85      0.796   4.08e+05      -82.3   4.83e-05      0.528          1          1
-       32          0       7.58       7.79       1.59   6.19e+05      -82.6      0.226      0.504          1          1
-       64          1       7.73       3.94       3.25   6.37e+07      -82.9      0.214      0.624          1          1
-      128          4       7.88       4.27       6.37   1.17e+08      -87.9      0.352      0.608          1          1
-      256         14       7.64       4.44         13    2.3e+08      -87.4      0.402       0.61          1          1
-      512         24       7.95       5.09       25.4   4.65e+08      -87.1      0.377      0.579          1          1
- 1.02e+03         33       7.91       5.73       50.5   9.25e+08      -86.4      0.371       0.56          1          1
- 2.05e+03         79       7.89       5.74        101   1.86e+09      -86.9      0.432       0.56          1          1
-  4.1e+03        181       7.97       5.71        201   3.66e+09      -87.2       0.43      0.559          1          1
+        2          0       4.56        3.8      0.327   4.03e+06       -216          0       0.73          1          1
+        4          0       3.75       3.29      0.217   1.63e+05       -108    1.8e-64      0.773          1          1
+        8          0       5.61       6.73      0.414    2.3e+05      -83.1   0.000173      0.602          1          1
+       16          0       6.78       7.85      0.803   4.08e+05      -82.3   4.83e-05      0.528          1          1
+       32          0       7.58       7.79        1.6   6.19e+05      -82.6      0.226      0.504          1          1
+       64          1       7.73       3.94       3.29   6.37e+07      -82.9      0.214      0.624          1          1
+      128          4       7.88       4.27       6.39   1.17e+08      -87.9      0.352      0.608          1          1
+      256         14       7.64       4.44       12.7    2.3e+08      -87.4      0.402       0.61          1          1
+      512         24       7.95       5.09       25.6   4.65e+08      -87.1      0.377      0.579          1          1
+ 1.02e+03         33       7.91       5.73         51   9.25e+08      -86.4      0.371       0.56          1          1
+ 2.05e+03         79       7.89       5.74        102   1.86e+09      -86.9      0.432       0.56          1          1
+  4.1e+03        181       7.97       5.71        202   3.66e+09      -87.2       0.43      0.559          1          1
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Model Comparison

Let's now display the constraints provided by each data set in a single corner plot

# Create a corner plot / pair plot.
 using CairoMakie: Makie
 using PairPlots
@@ -319,4 +319,4 @@
         PairPlots.MarginDensity(),
         PairPlots.Scatter()
     )
-)
Example block output

We see that the constraints provided by the PMA, the astrometry, and the radial velocity data all individually overlap, and agree with the joint model constraint. This is means that none of the datasets are in tension with each other, which might suggest an issue with the data or with the modelling assumptions (e.g. single planet).

+)Example block output

We see that the constraints provided by the PMA, the astrometry, and the radial velocity data all individually overlap, and agree with the joint model constraint. This is means that none of the datasets are in tension with each other, which might suggest an issue with the data or with the modelling assumptions (e.g. single planet).

diff --git a/dev/cEri-plot-grid.png b/dev/cEri-plot-grid.png index 7ae1fba6..a3866158 100644 Binary files a/dev/cEri-plot-grid.png and b/dev/cEri-plot-grid.png differ diff --git a/dev/c_Eri_straight_line-plot-grid.png b/dev/c_Eri_straight_line-plot-grid.png index b802c3d9..ab37d15b 100644 Binary files a/dev/c_Eri_straight_line-plot-grid.png and b/dev/c_Eri_straight_line-plot-grid.png differ diff --git a/dev/chains/index.html b/dev/chains/index.html index 085d6df6..ab8c1530 100644 --- a/dev/chains/index.html +++ b/dev/chains/index.html @@ -2,4 +2,4 @@ Chains · Octofitter.jl

Chains

This page describes the format of the Monte Carlo chains created by Octofitter.jl

The output of the samplers in Octofitter is an MCMCChains.Chains object

A column will be present for each variable in your model, both defined in the Priors blocks or as Derived variables.

Variables defined for the System as a whole can be accessed directly. For example:

chain["M"]

This will return an array of μ values from the posterior. The format is a matrix of N-samples by N-walkers.

Variables defined for an individual Planet are grouped according to the standard MCMCChains format. For example:

chain["b_a"]

This returns an array of semi-major axis values (a) for the planet b sampled from the posterior.

Diagnostics

Printing the chains will display a number of useful summaries for each quantity, like the mean, 0.25, 0.5, and 0.75 quantiles, and convergence metrics. See MCMCChains documentation for more details.

Exporting Chains

There are two useful ways to export chains. One is with the JLD2 library which preserves all the information and structure of the chain (but is only easy to open again in Julia) and the other is converting them to tables.

As a table

You can convert your chains to any Tables.jl compatible table. TypedTables.Table is included with this package, but DataFrames.DataFrame works well too.

tbl = Table(chain)
 
 using DataFrames
-df = DataFrame(chain)

You can then use a wide variety of Tables.jl source or sink libraries to persist your data to a file or database. The easiest is probably Arrow.jl:

Arrow.write("mychain.csv", tbl)

Other useful formats could be CSV.jl or SQLite.jl. In these formats, the data can be archived and imported easily into other programs; however, there is not yet an automatic way to return the data into the MCMCChains.Chains format it originated in.

+df = DataFrame(chain)

You can then use a wide variety of Tables.jl source or sink libraries to persist your data to a file or database. The easiest is probably Arrow.jl:

Arrow.write("mychain.csv", tbl)

Other useful formats could be CSV.jl or SQLite.jl. In these formats, the data can be archived and imported easily into other programs; however, there is not yet an automatic way to return the data into the MCMCChains.Chains format it originated in.

diff --git a/dev/compat-orbitize/index.html b/dev/compat-orbitize/index.html index 5ba22cb3..0e8af8c7 100644 --- a/dev/compat-orbitize/index.html +++ b/dev/compat-orbitize/index.html @@ -1,2 +1,2 @@ -Orbitize! Compatibility · Octofitter.jl

Compatibility with Orbitize!

The Orbitize! python library is a popular package for fitting astrometric and radial velocity orbits.

Octofitter has support for loading and saving posteriors in HDF5 format–the same format used by Orbitize!. This is useful if you want to load an Orbitize! posterior into Octofitter for plotting or to compare results. Similarily, you can export Octofitter chains to use with Orbitize! analysis tools, including the popular whereistheplanet.com website for predicting planet locations from stored posteriors.

Warning

The Orbitize! import/export functionality only works with simple models with visual orbits and only one companion.

In addition, it is possible to load a orbit posterior and/or astrometry data directly from whereistheplanet.com by target name.

Loading an Orbitize! posterior

chain = Octofitter.loadhdf5("fname.h5")

Save a posterior in Orbitize! format

Octofitter.savehdf5("fname.h5", model, chain)

Loading an Orbitize! posterior saved to Whereistheplanet.com

chain = Octofitter.loadhdf5("51erib",)

Loading Astrometry Data saved to Whereistheplanet.com

astrom_like1, astro_like2 = Octofitter.Whereistheplanet_astrom("51erib"; object=1)

Two different astrometry likelihood objects are returned since orbitize supports both PA/sep and RA/DEC formats. Octofitter also supports both formats, but they must be placed into separate likelihood objects. Simply add both to the model to include all data.

+Orbitize! Compatibility · Octofitter.jl

Compatibility with Orbitize!

The Orbitize! python library is a popular package for fitting astrometric and radial velocity orbits.

Octofitter has support for loading and saving posteriors in HDF5 format–the same format used by Orbitize!. This is useful if you want to load an Orbitize! posterior into Octofitter for plotting or to compare results. Similarily, you can export Octofitter chains to use with Orbitize! analysis tools, including the popular whereistheplanet.com website for predicting planet locations from stored posteriors.

Warning

The Orbitize! import/export functionality only works with simple models with visual orbits and only one companion.

In addition, it is possible to load a orbit posterior and/or astrometry data directly from whereistheplanet.com by target name.

Loading an Orbitize! posterior

chain = Octofitter.loadhdf5("fname.h5")

Save a posterior in Orbitize! format

Octofitter.savehdf5("fname.h5", model, chain)

Loading an Orbitize! posterior saved to Whereistheplanet.com

chain = Octofitter.loadhdf5("51erib",)

Loading Astrometry Data saved to Whereistheplanet.com

astrom_like1, astro_like2 = Octofitter.Whereistheplanet_astrom("51erib"; object=1)

Two different astrometry likelihood objects are returned since orbitize supports both PA/sep and RA/DEC formats. Octofitter also supports both formats, but they must be placed into separate likelihood objects. Simply add both to the model to include all data.

diff --git a/dev/cross-validation/index.html b/dev/cross-validation/index.html index 635017da..9c93b9d3 100644 --- a/dev/cross-validation/index.html +++ b/dev/cross-validation/index.html @@ -33,4 +33,4 @@ scatter!(ax, result.pointwise(:p_eff)) -fig +fig diff --git a/dev/custom-likelihood/index.html b/dev/custom-likelihood/index.html index 9b534025..89a28ef1 100644 --- a/dev/custom-likelihood/index.html +++ b/dev/custom-likelihood/index.html @@ -47,4 +47,4 @@ astrometry_table = Table() return MyLikelihood(epoch=epochs, ra=ras, dec=decs) -end +end diff --git a/dev/derived/index.html b/dev/derived/index.html index 9478b760..fbd95ff6 100644 --- a/dev/derived/index.html +++ b/dev/derived/index.html @@ -45,4 +45,4 @@ M ~ Normal(45., 0.02) i ~ Normal(0.1, deg2rad(30.)) Ω ~ Normal(0.0, deg2rad(30.)) -end b c

Notice how i and Ω are defined as variables on the System. The two planets B & C instead just take their values from the System. This way we can enforce co-planarity between planets without e.g. rejection sampling.

Resolution Order

The order that variables are resolved is as follows:

You can use one derived variable from another based on their order in the @system or @planet block. You cannot access variables from a different planet inside a @planet block. If you need to do this, move the variable up to the @system block.

+end b c

Notice how i and Ω are defined as variables on the System. The two planets B & C instead just take their values from the System. This way we can enforce co-planarity between planets without e.g. rejection sampling.

Resolution Order

The order that variables are resolved is as follows:

You can use one derived variable from another based on their order in the @system or @planet block. You cannot access variables from a different planet inside a @planet block. If you need to do this, move the variable up to the @system block.

diff --git a/dev/extract-phot-astrom/index.html b/dev/extract-phot-astrom/index.html index 966e088f..b3ba678e 100644 --- a/dev/extract-phot-astrom/index.html +++ b/dev/extract-phot-astrom/index.html @@ -50,8 +50,8 @@ Iterations = 1:1:10000 Number of chains = 1 Samples per chain = 10000 -Wall duration = 2.76 seconds -Compute duration = 2.76 seconds +Wall duration = 2.98 seconds +Compute duration = 2.98 seconds parameters = plx, b_sep, b_pa, b_L internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -60,9 +60,9 @@ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ plx 24.4620 0.0000 0.0000 NaN NaN NaN ⋯ - b_sep 1704.9364 5.4858 0.6542 195.7072 52.6697 1.0006 ⋯ - b_pa 1.2330 0.0020 0.0000 3662.8178 4922.1008 1.0005 ⋯ - b_L 0.0001 0.0000 0.0000 2449.4307 3502.3575 1.0002 ⋯ + b_sep 1705.1620 5.7527 0.6542 194.1852 74.1903 1.0027 ⋯ + b_pa 1.2330 0.0020 0.0000 5334.6849 6388.9112 0.9999 ⋯ + b_L 0.0001 0.0000 0.0000 2947.3812 3405.4661 1.0000 ⋯ 1 column omitted Quantiles @@ -70,12 +70,12 @@ Symbol Float64 Float64 Float64 Float64 Float64 plx 24.4620 24.4620 24.4620 24.4620 24.4620 - b_sep 1698.5460 1701.8755 1703.6383 1706.3640 1723.5214 - b_pa 1.2296 1.2315 1.2328 1.2344 1.2370 + b_sep 1698.6201 1701.9089 1703.7075 1706.5966 1723.8912 + b_pa 1.2297 1.2314 1.2328 1.2345 1.2370 b_L 0.0001 0.0001 0.0001 0.0001 0.0001

Sample from the model (globally)

You could also try sampling across the entire image, without necessarily specifying a starting position. Note that if there are multiple candidates, taking the naive mean and standard deviation will average across all planets.

using Pigeons
 model.starting_points = nothing # reset starting points
 chain, pt = octofit_pigeons(model, n_rounds=11)
(chain = MCMC chain (2048×8×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

Access results

samples_sep = chain[:b_sep]
 samples_pa = chain[:b_pa]
 println("The median separation is ", median(samples_sep))
The median separation is 1703.7929269387687

Visualize

using CairoMakie, PairPlots
-octocorner(model,chain)
Example block output +octocorner(model,chain)Example block output diff --git a/dev/fit-coplanar/b21f8340.png b/dev/fit-coplanar/b21f8340.png deleted file mode 100644 index 8a846134..00000000 Binary files a/dev/fit-coplanar/b21f8340.png and /dev/null differ diff --git a/dev/fit-coplanar/dcc95f9b.png b/dev/fit-coplanar/dcc95f9b.png new file mode 100644 index 00000000..57039a7b Binary files /dev/null and b/dev/fit-coplanar/dcc95f9b.png differ diff --git a/dev/fit-coplanar/index.html b/dev/fit-coplanar/index.html index e43cf4ca..1119e525 100644 --- a/dev/fit-coplanar/index.html +++ b/dev/fit-coplanar/index.html @@ -99,20 +99,20 @@ ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ) ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── - 2 0 3.54 3.03 0.455 6.97e+06 -1.83e+07 0 0.788 1 1 - 4 0 4.62 5.1 0.388 1.57e+05 -1.06e+06 0 0.686 1 1 - 8 0 4.97 5.15 0.746 2.64e+05 -6.87e+05 0 0.674 1 1 - 16 0 6.47 6.08 1.43 4.11e+05 -8.53e+04 0 0.595 1 1 - 32 0 7.21 7.26 2.85 6.38e+05 -1.2e+04 0 0.533 1 1 - 64 0 7.46 5.3 5.7 7.1e+07 -214 0 0.588 1 1 - 128 3 8.5 5.75 11.2 1.3e+08 -232 1.16e-185 0.54 1 1 - 256 0 9.3 6.8 21.5 2.56e+08 -230 0 0.481 1 1 - 512 4 9.99 8.32 41.9 5.11e+08 -207 8.98e-20 0.409 1 1 - 1.02e+03 7 9.86 8.3 79.4 9.78e+08 -203 4.6e-26 0.414 1 1 -────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Plots the orbits:

octoplot(model, results)
Example block output

Corner plot:

octocorner(model, results, small=true)
Example block output

Now examine the period ratio:

hist(
+        2          0       3.54       3.03       0.46   6.97e+06  -1.83e+07          0      0.788          1          1
+        4          0       4.62        5.1      0.389   1.57e+05  -1.06e+06          0      0.686          1          1
+        8          0       4.97       5.15      0.748   2.64e+05  -6.87e+05          0      0.674          1          1
+       16          0       6.47       6.08       1.44   4.11e+05  -8.53e+04          0      0.595          1          1
+       32          0       7.21       7.26       2.84   6.38e+05   -1.2e+04          0      0.533          1          1
+       64          0       7.46        5.3       5.61    7.1e+07       -214          0      0.588          1          1
+      128          3        8.5       5.75         11    1.3e+08       -232  1.16e-185       0.54          1          1
+      256          0        9.3        6.8       21.3   2.56e+08       -230          0      0.481          1          1
+      512          4       9.99       8.32       41.2   5.11e+08       -207   8.98e-20      0.409          1          1
+ 1.02e+03          7       9.86        8.3       78.2   9.78e+08       -203    4.6e-26      0.414          1          1
+────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Plots the orbits:

octoplot(model, results)
Example block output

Corner plot:

octocorner(model, results, small=true)
Example block output

Now examine the period ratio:

hist(
     results[:b_P][:] ./ results[:c_P][:],
     axis=(;
         xlabel="period ratio",
         ylabel="counts",
     )
-)
Example block output +)Example block output diff --git a/dev/fit-grav-wide/index.html b/dev/fit-grav-wide/index.html index 4054bdfc..264286e9 100644 --- a/dev/fit-grav-wide/index.html +++ b/dev/fit-grav-wide/index.html @@ -104,4 +104,4 @@ K ~ Uniform(0, 1) # contrast ratio variable K_band_spectrum = fill(b.K, $(length(vis_like.table.eff_wave[1]))) -end +end diff --git a/dev/fit-interfere/03d5627d.png b/dev/fit-interfere/03d5627d.png new file mode 100644 index 00000000..52cb8495 Binary files /dev/null and b/dev/fit-interfere/03d5627d.png differ diff --git a/dev/fit-interfere/0c5e6bb5.png b/dev/fit-interfere/0c5e6bb5.png deleted file mode 100644 index 05464ac7..00000000 Binary files a/dev/fit-interfere/0c5e6bb5.png and /dev/null differ diff --git a/dev/fit-interfere/index.html b/dev/fit-interfere/index.html index d4fa59c2..495d1971 100644 --- a/dev/fit-interfere/index.html +++ b/dev/fit-interfere/index.html @@ -3,7 +3,7 @@ using OctofitterInterferometry using Distributions using CairoMakie -using PairPlots
┌ Warning: Module Octofitter with build ID fafbfcfd-c563-0173-0000-003d53fb7590 is missing from the cache.
+using PairPlots
┌ Warning: Module Octofitter with build ID fafbfcfd-1710-b70a-0000-003b9bab10ca is missing from the cache.
 │ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
 └ @ Base loading.jl:1948

Download simulated JWST AMI observations from our examples folder on GitHub:

download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2023_1_.oifits", "Sim_data_2023_1_.oifits")
 download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2023_2_.oifits", "Sim_data_2023_2_.oifits")
@@ -96,19 +96,19 @@
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
 ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
-        2          0        5.2       3.53      0.424   8.63e+07  -6.22e+04          0      0.718          1          1
-        4          0       3.88        5.2      0.443   1.58e+08        158          0      0.707          1          1
-        8          0       4.53       6.09       0.84   3.07e+08        176    3.1e-05      0.658          1          1
-       16          0        6.1       5.57       1.66   6.07e+08        178   2.87e-24      0.624          1          1
-       32          0        5.9       5.86       3.32   1.22e+09        178    0.00425      0.621          1          1
-       64          6       6.49       2.71       6.94   2.56e+09        176      0.128      0.703          1          1
+        2          0        5.2       3.53      0.405   8.63e+07  -6.22e+04          0      0.718          1          1
+        4          0       3.88        5.2       0.43   1.58e+08        158          0      0.707          1          1
+        8          0       4.53       6.09      0.851   3.07e+08        176    3.1e-05      0.658          1          1
+       16          0        6.1       5.57       1.63   6.07e+08        178   2.87e-24      0.624          1          1
+       32          0        5.9       5.86       3.31   1.22e+09        178    0.00425      0.621          1          1
+       64          6       6.49       2.71       6.93   2.56e+09        176      0.128      0.703          1          1
       128         15       6.42       2.75       13.5   5.13e+09        176      0.188      0.704          1          1
       256         29       6.78       2.61       26.9   1.02e+10        176      0.344      0.697          1          1
       512         71       6.63       2.68       53.8   2.05e+10        176      0.481        0.7          1          1
  1.02e+03        135       6.82        2.7        108    4.1e+10        176      0.507      0.693          1          1
 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Note that we use Pigeons paralell tempered sampling (octofit_pigeons) instead of HMC (octofit) because interferometry data is almost always multi-modal (or more precisely non-convex, there is often still a single mode that dominates).

Examine the recovered photometry posterior:

hist(results[:b_contrast_F480M][:], axis=(;xlabel="F480M"))
Example block output

Determine the significance of the detection:

using Statistics
 phot = results[:b_contrast_F480M][:]
-snr = mean(phot)/std(phot)
7.290535815433398

Plot the resulting orbit:

octoplot(model, results)
Example block output

Plot only the position at each epoch:

using PlanetOrbits
+snr = mean(phot)/std(phot)
7.290535815433398

Plot the resulting orbit:

octoplot(model, results)
Example block output

Plot only the position at each epoch:

using PlanetOrbits
 els = Octofitter.construct_elements(results,:b,:);
 fig = Makie.Figure()
 ax = Makie.Axis(
@@ -130,4 +130,4 @@
 Makie.Legend(fig[1,2], ax, "date")
 fig
Example block output

Finally we can examine the joint photometry and orbit posterior as a corner plot:

using PairPlots
 using CairoMakie: Makie
-octocorner(model, results)
Example block output +octocorner(model, results)Example block output diff --git a/dev/fit-likemap/7b07afd2.png b/dev/fit-likemap/7b07afd2.png deleted file mode 100644 index 8b36cc3f..00000000 Binary files a/dev/fit-likemap/7b07afd2.png and /dev/null differ diff --git a/dev/fit-likemap/ee4296e0.png b/dev/fit-likemap/ee4296e0.png new file mode 100644 index 00000000..c1fd39df Binary files /dev/null and b/dev/fit-likemap/ee4296e0.png differ diff --git a/dev/fit-likemap/index.html b/dev/fit-likemap/index.html index d07f938c..0fa8e0b3 100644 --- a/dev/fit-likemap/index.html +++ b/dev/fit-likemap/index.html @@ -144,17 +144,17 @@ ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ) ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── - 2 0 2.38 3.52 0.255 6.97e+06 -20.8 1.21e-25 0.81 1 1 - 4 0 3.85 4.45 0.0596 1.53e+05 -23.9 9.31e-07 0.732 1 1 - 8 0 4.19 3.74 0.118 2.7e+05 -22.4 0.2 0.744 1 1 - 16 0 4.77 4.64 0.234 4.54e+05 -21.5 0.221 0.696 1 1 - 32 0 4.23 4.83 0.467 7.27e+05 -21 0.471 0.708 1 1 - 64 4 4.93 3.08 1.07 4.65e+07 -20.6 0.536 0.742 1 1 - 128 16 4.88 3.08 1.82 7.93e+07 -21 0.511 0.743 1 1 - 256 36 4.71 3.09 3.69 1.58e+08 -21.1 0.623 0.748 1 1 - 512 77 4.75 3.09 7.39 3.16e+08 -20.9 0.634 0.747 1 1 - 1.02e+03 164 4.82 3.09 14.7 6.28e+08 -21 0.61 0.745 1 1 -────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Note

octofit_pigeons scales very well across multiple cores. Start julia with julia --threads=auto to make sure you have multiple threads available for sampling.

Display the results:

octoplot(model, chain)
Example block output

Corner plot:

using CairoMakie, PairPlots
+        2          0       2.38       3.52      0.283   6.97e+06      -20.8   1.21e-25       0.81          1          1
+        4          0       3.85       4.45     0.0601   1.53e+05      -23.9   9.31e-07      0.732          1          1
+        8          0       4.19       3.74      0.119    2.7e+05      -22.4        0.2      0.744          1          1
+       16          0       4.77       4.64      0.235   4.54e+05      -21.5      0.221      0.696          1          1
+       32          0       4.23       4.83       0.47   7.27e+05        -21      0.471      0.708          1          1
+       64          4       4.93       3.08        1.1   4.65e+07      -20.6      0.536      0.742          1          1
+      128         16       4.88       3.08       1.86   7.93e+07        -21      0.511      0.743          1          1
+      256         36       4.71       3.09       3.74   1.58e+08      -21.1      0.623      0.748          1          1
+      512         77       4.75       3.09       7.47   3.16e+08      -20.9      0.634      0.747          1          1
+ 1.02e+03        164       4.82       3.09       14.9   6.28e+08        -21       0.61      0.745          1          1
+────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Note

octofit_pigeons scales very well across multiple cores. Start julia with julia --threads=auto to make sure you have multiple threads available for sampling.

Display the results:

octoplot(model, chain)
Example block output

Corner plot:

using CairoMakie, PairPlots
 octocorner(model,chain,small=true,)
Example block output

And finally let's look at the posterior predictive distributions at both epochs:

els = Octofitter.construct_elements(chain,:b, :)
 x = raoff.(els, loglikemap.table.epoch[1])
 y = decoff.(els, loglikemap.table.epoch[1])
@@ -182,4 +182,4 @@
     )
 )
Example block output

Resume sampling for additional rounds

If you would like to add additional rounds of sampling, you may do the following:

pt = increment_n_rounds!(pt, 2)
 chain, pt = octofit_pigeons(pt)
(chain = MCMC chain (4096×19×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

Updated corner plot:

using CairoMakie, PairPlots
-octocorner(model,chain,small=false,)
Example block output +octocorner(model,chain,small=false,)Example block output diff --git a/dev/fit-rv-astrom/index.html b/dev/fit-rv-astrom/index.html index d81ac20d..53f36472 100644 --- a/dev/fit-rv-astrom/index.html +++ b/dev/fit-rv-astrom/index.html @@ -73,8 +73,8 @@ Iterations = 1:1:5000 Number of chains = 1 Samples per chain = 5000 -Wall duration = 18.86 seconds -Compute duration = 18.86 seconds +Wall duration = 19.21 seconds +Compute duration = 19.21 seconds parameters = M, jitter1, plx, b_e, b_a, b_mass, b_i, b_Ωy, b_Ωx, b_ωy, b_ωx, b_θy, b_θx, b_Ω, b_ω, b_θ, b_tp internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -122,4 +122,4 @@ b_ω -3.0731 -2.1397 0.2804 2.5150 3.0533 b_θ -2.9861 -2.9788 -2.9748 -2.9707 -2.9640 b_tp 57670.8627 58368.7893 58659.7674 58752.8917 58799.7988 -

Display results as a corner plot:

octocorner(model,results, small=true)
Example block output

Plot RV curve, phase folded curve, and binned residuals:

Octofitter.rvpostplot(model, results,)
Example block output

Display RV, PMA, astrometry, relative separation, position angle, and 3D projected views:

octoplot(model, results)
Example block output +

Display results as a corner plot:

octocorner(model,results, small=true)
Example block output

Plot RV curve, phase folded curve, and binned residuals:

Octofitter.rvpostplot(model, results,)
Example block output

Display RV, PMA, astrometry, relative separation, position angle, and 3D projected views:

octoplot(model, results)
Example block output diff --git a/dev/fit-rv-rel/index.html b/dev/fit-rv-rel/index.html index 396bda5d..06d72184 100644 --- a/dev/fit-rv-rel/index.html +++ b/dev/fit-rv-rel/index.html @@ -78,8 +78,8 @@ Iterations = 1:1:1000 Number of chains = 1 Samples per chain = 1000 -Wall duration = 2.82 seconds -Compute duration = 2.82 seconds +Wall duration = 2.88 seconds +Compute duration = 2.88 seconds parameters = M, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_τy, b_τx, b_gamma, b_ω, b_Ω, b_τ, b_P, b_tp internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -125,4 +125,4 @@ b_τ -0.4725 -0.2589 -0.0553 0.2620 0.4771 b_P 1.2774 1.3042 1.3204 1.3365 1.3692 b_tp 5769.0376 5875.6452 5973.3538 6126.1830 6229.2390 -
octoplot(model, chain, ts=range(4500, 8000, length=200), show_physical_orbit=true)
Example block output +
octoplot(model, chain, ts=range(4500, 8000, length=200), show_physical_orbit=true)
Example block output diff --git a/dev/hipparcos/7ed44b92.png b/dev/hipparcos/7ed44b92.png deleted file mode 100644 index b802c3d9..00000000 Binary files a/dev/hipparcos/7ed44b92.png and /dev/null differ diff --git a/dev/hipparcos/870603d6.png b/dev/hipparcos/870603d6.png deleted file mode 100644 index 50704d21..00000000 Binary files a/dev/hipparcos/870603d6.png and /dev/null differ diff --git a/dev/hipparcos/b39a95b2.png b/dev/hipparcos/b39a95b2.png new file mode 100644 index 00000000..ab37d15b Binary files /dev/null and b/dev/hipparcos/b39a95b2.png differ diff --git a/dev/hipparcos/bb5df9ea.png b/dev/hipparcos/bb5df9ea.png new file mode 100644 index 00000000..8bf308fe Binary files /dev/null and b/dev/hipparcos/bb5df9ea.png differ diff --git a/dev/hipparcos/dcb2f900.png b/dev/hipparcos/dcb2f900.png deleted file mode 100644 index 7ae1fba6..00000000 Binary files a/dev/hipparcos/dcb2f900.png and /dev/null differ diff --git a/dev/hipparcos/eedede49.png b/dev/hipparcos/eedede49.png new file mode 100644 index 00000000..a3866158 Binary files /dev/null and b/dev/hipparcos/eedede49.png differ diff --git a/dev/hipparcos/index.html b/dev/hipparcos/index.html index 9fc1cdb8..c47d825a 100644 --- a/dev/hipparcos/index.html +++ b/dev/hipparcos/index.html @@ -42,8 +42,8 @@ Iterations = 1:1:4000 Number of chains = 1 Samples per chain = 4000 -Wall duration = 1.19 seconds -Compute duration = 1.19 seconds +Wall duration = 1.17 seconds +Compute duration = 1.17 seconds parameters = plx, pmra, pmdec, ra_hip_offset_mas, dec_hip_offset_mas, M, rv, dec, ra, ref_epoch, b_mass, b_e, b_ω, b_a, b_i, b_Ω, b_tp internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -51,15 +51,15 @@ parameters mean std mcse ess_bulk ess_tail ⋯ Symbol Float64 Float64 Float64 Float64 Float64 ⋯ - plx 95.8653 0.0000 NaN NaN NaN ⋯ - pmra -55.7650 0.0000 0.0000 NaN NaN ⋯ - pmdec -48.6715 0.0000 0.0000 NaN NaN ⋯ - ra_hip_offset_mas 393.7259 0.0000 0.0000 NaN NaN ⋯ - dec_hip_offset_mas -1293.3630 0.0000 0.0000 NaN NaN ⋯ + plx 71.3895 0.0000 0.0000 NaN NaN ⋯ + pmra 60.9473 0.0000 0.0000 NaN NaN ⋯ + pmdec 23.8452 0.0000 0.0000 NaN NaN ⋯ + ra_hip_offset_mas 158.3912 0.0000 0.0000 NaN NaN ⋯ + dec_hip_offset_mas 1169.3996 0.0000 0.0000 NaN NaN ⋯ M 1.0000 0.0000 NaN NaN NaN ⋯ rv 0.0000 0.0000 NaN NaN NaN ⋯ dec -2.4733 0.0000 0.0000 NaN NaN ⋯ - ra 69.4001 0.0000 0.0000 NaN NaN ⋯ + ra 69.4008 0.0000 0.0000 NaN NaN ⋯ ref_epoch 48348.5625 0.0000 NaN NaN NaN ⋯ b_mass 0.0000 0.0000 NaN NaN NaN ⋯ b_e 0.0000 0.0000 NaN NaN NaN ⋯ @@ -74,15 +74,15 @@ parameters 2.5% 25.0% 50.0% 75.0% ⋯ Symbol Float64 Float64 Float64 Float64 ⋯ - plx 95.8653 95.8653 95.8653 95.8653 ⋯ - pmra -55.7650 -55.7650 -55.7650 -55.7650 - ⋯ - pmdec -48.6715 -48.6715 -48.6715 -48.6715 - ⋯ - ra_hip_offset_mas 393.7259 393.7259 393.7259 393.7259 3 ⋯ - dec_hip_offset_mas -1293.3630 -1293.3630 -1293.3630 -1293.3630 -12 ⋯ + plx 71.3895 71.3895 71.3895 71.3895 ⋯ + pmra 60.9473 60.9473 60.9473 60.9473 ⋯ + pmdec 23.8452 23.8452 23.8452 23.8452 ⋯ + ra_hip_offset_mas 158.3912 158.3912 158.3912 158.3912 1 ⋯ + dec_hip_offset_mas 1169.3996 1169.3996 1169.3996 1169.3996 11 ⋯ M 1.0000 1.0000 1.0000 1.0000 ⋯ rv 0.0000 0.0000 0.0000 0.0000 ⋯ dec -2.4733 -2.4733 -2.4733 -2.4733 ⋯ - ra 69.4001 69.4001 69.4001 69.4001 ⋯ + ra 69.4008 69.4008 69.4008 69.4008 ⋯ ref_epoch 48348.5625 48348.5625 48348.5625 48348.5625 483 ⋯ b_mass 0.0000 0.0000 0.0000 0.0000 ⋯ b_e 0.0000 0.0000 0.0000 0.0000 ⋯ @@ -92,7 +92,7 @@ b_Ω 0.0000 0.0000 0.0000 0.0000 ⋯ b_tp 0.0000 0.0000 0.0000 0.0000 ⋯ 1 column omitted -

Plot the posterior values:

octoplot(model,chain,show_astrom=false,show_astrom_time=false)
Example block output

We now visualize the model fit compared to the Hipparcos catalog values:

using LinearAlgebra, StatsBase
+

Plot the posterior values:

octoplot(model,chain,show_astrom=false,show_astrom_time=false)
Example block output

We now visualize the model fit compared to the Hipparcos catalog values:

using LinearAlgebra, StatsBase
 fig = Figure(size=(1080,720))
 j = i = 1
 for prop in (
@@ -133,7 +133,7 @@
     lines!(ax, nxs, pdf.(n,nxs), label="Hipparcos Catalog", color=:black, linewidth=2)
 end
 Legend(fig[i-1,j+1],ax,tellwidth=false)
-fig
Example block output

Constrain Planet Mass

We now allow the planet to have a non zero mass and have free orbit. We start by retrieving relative astrometry data on the planet, collated by Jason Wang and co. on whereistheplanet.com.

astrom_like1,astrom_like2 = Octofitter.Whereistheplanet_astrom("51erib",object=1)

We specify our full model:

@planet b AbsoluteVisual{KepOrbit} begin
+fig
Example block output

Constrain Planet Mass

We now allow the planet to have a non zero mass and have free orbit. We start by retrieving relative astrometry data on the planet, collated by Jason Wang and co. on whereistheplanet.com.

astrom_like1,astrom_like2 = Octofitter.Whereistheplanet_astrom("51erib",object=1)

We specify our full model:

@planet b AbsoluteVisual{KepOrbit} begin
     a ~ truncated(Normal(10,1),lower=0.1)
     e ~ Uniform(0,0.99)
     ω ~ UniformCircular()
@@ -171,8 +171,8 @@
 Iterations        = 1:1:256
 Number of chains  = 1
 Samples per chain = 256
-Wall duration     = 439.51 seconds
-Compute duration  = 439.51 seconds
+Wall duration     = 440.66 seconds
+Compute duration  = 440.66 seconds
 parameters        = M_pri, M_sec, plx, pmra, pmdec, ra_hip_offset_mas, dec_hip_offset_mas, M, rv, dec, ra, ref_epoch, b_a, b_e, b_ωy, b_ωx, b_i, b_Ωy, b_Ωx, b_θy, b_θx, b_ω, b_Ω, b_θ, b_tp, b_mass
 internals         = loglike, logpost, logprior, pigeons_logpotential
 
@@ -223,4 +223,4 @@
                  b_i       2.3387       2.3997       2.4350       2.4660       ⋯
           ⋮                ⋮            ⋮            ⋮            ⋮            ⋱
                                                      1 column and 9 rows omitted
-
octoplot(model, chain, show_mass=true)
Example block output

We see that we constrained both the orbit and the parallax. The mass is not strongly constrained by Hipparcos.

+
octoplot(model, chain, show_mass=true)
Example block output

We see that we constrained both the orbit and the parallax. The mass is not strongly constrained by Hipparcos.

diff --git a/dev/images/17dabbdd.png b/dev/images/17dabbdd.png new file mode 100644 index 00000000..5153a77a Binary files /dev/null and b/dev/images/17dabbdd.png differ diff --git a/dev/images/31919b12.png b/dev/images/31919b12.png new file mode 100644 index 00000000..327626e8 Binary files /dev/null and b/dev/images/31919b12.png differ diff --git a/dev/images/48c3113e.png b/dev/images/48c3113e.png new file mode 100644 index 00000000..b5c62add Binary files /dev/null and b/dev/images/48c3113e.png differ diff --git a/dev/images/7d82e334.png b/dev/images/7d82e334.png deleted file mode 100644 index f65dcc3a..00000000 Binary files a/dev/images/7d82e334.png and /dev/null differ diff --git a/dev/images/c8eb09a7.png b/dev/images/c8eb09a7.png deleted file mode 100644 index 31b5a8f5..00000000 Binary files a/dev/images/c8eb09a7.png and /dev/null differ diff --git a/dev/images/ce7df5c8.png b/dev/images/ce7df5c8.png deleted file mode 100644 index 2f31ac27..00000000 Binary files a/dev/images/ce7df5c8.png and /dev/null differ diff --git a/dev/images/index.html b/dev/images/index.html index b4051539..947ad791 100644 --- a/dev/images/index.html +++ b/dev/images/index.html @@ -123,16 +123,16 @@ ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ) ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── - 2 0 5.69 3.78 0.252 4.03e+06 33.7 2.07e-09 0.694 1 1 - 4 0 2.03 3.58 0.126 1.61e+05 58.8 0.00817 0.819 1 1 - 8 0 4.46 4.83 0.232 2.59e+05 64.7 0.0682 0.7 1 1 - 16 0 4.54 5.53 0.485 4.34e+05 61.9 1.54e-07 0.675 1 1 - 32 0 6.14 5.35 0.965 6.64e+05 60.7 0.191 0.629 1 1 - 64 3 5.3 3.16 2.3 5.73e+07 64.4 0.248 0.727 1 1 - 128 2 5.15 4.16 4.06 1e+08 64.6 0.26 0.7 1 1 - 256 14 4.96 4.59 8.18 1.95e+08 65.8 0.308 0.692 1 1 - 512 25 5.2 4.7 16.2 3.92e+08 65.1 0.472 0.681 1 1 - 1.02e+03 63 5.28 4.7 31.9 7.69e+08 65.4 0.551 0.678 1 1 + 2 0 5.69 3.78 0.271 4.03e+06 33.7 2.07e-09 0.694 1 1 + 4 0 2.03 3.58 0.127 1.61e+05 58.8 0.00817 0.819 1 1 + 8 0 4.46 4.83 0.235 2.59e+05 64.7 0.0682 0.7 1 1 + 16 0 4.54 5.53 0.483 4.34e+05 61.9 1.54e-07 0.675 1 1 + 32 0 6.14 5.35 0.975 6.64e+05 60.7 0.191 0.629 1 1 + 64 3 5.3 3.16 2.24 5.73e+07 64.4 0.248 0.727 1 1 + 128 2 5.15 4.16 4.11 1e+08 64.6 0.26 0.7 1 1 + 256 14 4.96 4.59 8.22 1.95e+08 65.8 0.308 0.692 1 1 + 512 25 5.2 4.7 16.4 3.92e+08 65.1 0.472 0.681 1 1 + 1.02e+03 63 5.28 4.7 32.3 7.69e+08 65.4 0.551 0.678 1 1 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Note

octofit_pigeons scales very well across multiple cores. Start julia with julia --threads=auto to make sure you have multiple threads available for sampling.

Diagnostics

The first thing you should do with your results is check a few diagnostics to make sure the sampler converged as intended.

The acceptance rate should be somewhat lower than when fitting just astrometry, e.g. around the 0.6 target.

You can make a trace plot:

lines(
     chain["b_a"][:],
     axis=(;
@@ -146,7 +146,7 @@
         xlabel="lag",
         ylabel="autocorrelation",
     )
-)
Example block output

For this model, there is somewhat higher correlation between samples. Some thinning to remove this correlation is recommended.

Analysis

We can now view the orbit fit:

fig = octoplot(model, chain)
Example block output

With a bit of work, we can plot one of the images under the orbit.

fig = octoplot(model, chain)
+)
Example block output

For this model, there is somewhat higher correlation between samples. Some thinning to remove this correlation is recommended.

Analysis

We can now view the orbit fit:

fig = octoplot(model, chain)
Example block output

With a bit of work, we can plot one of the images under the orbit.

fig = octoplot(model, chain)
 ax = fig.content[1] # grap first axis in the figure
 
 # We have to do some annoying work to get the image orientated correctly
@@ -163,7 +163,7 @@
 Colorbar(fig[1,2], h, label="image flux")
 
 Makie.resize_to_layout!(fig)
-fig
Example block output

Another useful view would be the orbits over a stack of the maximum pixel values of all images.

fig = octoplot(model, chain)
+fig
Example block output

Another useful view would be the orbits over a stack of the maximum pixel values of all images.

fig = octoplot(model, chain)
 ax = fig.content[1] # grap first axis in the figure
 
 # We have to do some annoying work to get the image orientated correctly
@@ -176,6 +176,6 @@
 h = heatmap!(ax, imgax1, imgax2, collect(img), colormap=:greys)
 Makie.translate!(h, 0,0,-1) # Send heatmap to back of the plot
 Makie.resize_to_layout!(fig)
-fig
Example block output

Pair Plot

We can show the relationships between variables on a pair plot (aka corner plot):

using CairoMakie, PairPlots
+fig
Example block output

Pair Plot

We can show the relationships between variables on a pair plot (aka corner plot):

using CairoMakie, PairPlots
 octocorner(model, chain, small=true)
Example block output

Note that this time, we also show the recovered photometry in the corner plot.

Assessing Detections

To assess a detection, we can treat all the orbital variables as nuisance parameters. We start by plotting the marginal distribution of the flux parameter, H:

hist(chain["b_H"][:], axis=(xlabel="H", ylabel="counts"))
Example block output

We can calculate an analog of the traditional signal to noise ratio (SNR) using that same histogram:

flux = chain["b_H"]
-snr = mean(flux)/std(flux)
15.172535840305645

It might be better to consider a related measure, like the median flux over the interquartile distance. This will depend on your application.

+snr = mean(flux)/std(flux)
15.172535840305645

It might be better to consider a related measure, like the median flux over the interquartile distance. This will depend on your application.

diff --git a/dev/index.html b/dev/index.html index a6729b44..bfeefbbe 100644 --- a/dev/index.html +++ b/dev/index.html @@ -16,4 +16,4 @@ title = {{PairPlots.jl} Beautiful and flexible visualizations of high dimensional data}, year = {2023}, howpublished = {\url{https://sefffal.github.io/PairPlots.jl/dev}}, -}

Ready?

Ready to get started? Follow our installation guide and then follow our first tutorial.

+}

Ready?

Ready to get started? Follow our installation guide and then follow our first tutorial.

diff --git a/dev/installation/index.html b/dev/installation/index.html index 0212bb59..2fe2773c 100644 --- a/dev/installation/index.html +++ b/dev/installation/index.html @@ -1,4 +1,4 @@ Installation · Octofitter.jl

Installation

The first step to using Octofitter.jl is to install Julia. If you're used to Python, don't worry –- Julia is easy to install, and you won't need to code anything other than changing your input data.

Installing Julia

Visit the julialang.org Downloads page, and select the latest stable version for your operating system. This is 1.10.0 at the time of writing. Click the [help] links next to your operating system if you require more detailed instructions.

Installing Octofitter

  1. Start julia in a terminal by running julia
  2. Type ] to enter package-mode (see Julia documentation for more details)
  3. Type add Octofitter Distributions CairoMakie PairPlots

You will need the Distributions,jl package so that you can specify priors for different parameters in your models. CairoMakie.jl is used for generating plots and isn't needed if you only want text-based summary outputs. PairPlots.jl (in combination with CairoMakie) is used for generating corner plots and can also be skipped if these aren't of interest.

Extension Packages

Some Octofitter functionality exists in extension packages, including radial velocity fitting. If you need one of these packages you can install them like so:

pkg> add OctofitterRadialVelocity
 pkg> add http://github.com/sefffal/Octofitter.jl:OctofitterImages
-pkg> add http://github.com/sefffal/Octofitter.jl:OctofitterInterferometry

These aren't included by default since they may include a number of heavier dependencies that aren't needed by all users. They are descibed further in relevant sections of the documentation.

Fitting your first model

Start with the Quick Start tutorial. It shows how one can model the orbit of one planet based on relative astrometry points.

+pkg> add http://github.com/sefffal/Octofitter.jl:OctofitterInterferometry

These aren't included by default since they may include a number of heavier dependencies that aren't needed by all users. They are descibed further in relevant sections of the documentation.

Fitting your first model

Start with the Quick Start tutorial. It shows how one can model the orbit of one planet based on relative astrometry points.

diff --git a/dev/kepler/index.html b/dev/kepler/index.html index 8ac97d55..cedb9655 100644 --- a/dev/kepler/index.html +++ b/dev/kepler/index.html @@ -1,2 +1,2 @@ -Kepler Solver · Octofitter.jl

Kepler Solver

The heart of this package is being able to take a set of Keplerian elements and output relative positions, velocities, etc. For this, we use PlanetOrbits.jl which adopts the same conventions as Orbitize!.

The Kepler solver used to go from mean anomaly to eccentric anomaly is a tweaked version copied from AstroLib.jl.

From AstroLib.jl:

Many different numerical methods exist to solve Kepler's equation. This function implements the algorithm proposed in Markley (1995) Celestial Mechanics and Dynamical Astronomy, 63, 101 (DOI:10.1007/BF00691917). This method is not iterative, requires only four transcendental function evaluations, and has been proved to be fast and efficient over the entire range of elliptic motion 0≤e≤10.

On my laptop, this solves for a single eccentric anomaly in just 47 ns. Since it is implemented in pure Julia, there is no overhead from calling into a C or Cython compiled function and no need for vectorization.

+Kepler Solver · Octofitter.jl

Kepler Solver

The heart of this package is being able to take a set of Keplerian elements and output relative positions, velocities, etc. For this, we use PlanetOrbits.jl which adopts the same conventions as Orbitize!.

The Kepler solver used to go from mean anomaly to eccentric anomaly is a tweaked version copied from AstroLib.jl.

From AstroLib.jl:

Many different numerical methods exist to solve Kepler's equation. This function implements the algorithm proposed in Markley (1995) Celestial Mechanics and Dynamical Astronomy, 63, 101 (DOI:10.1007/BF00691917). This method is not iterative, requires only four transcendental function evaluations, and has been proved to be fast and efficient over the entire range of elliptic motion 0≤e≤10.

On my laptop, this solves for a single eccentric anomaly in just 47 ns. Since it is implemented in pure Julia, there is no overhead from calling into a C or Cython compiled function and no need for vectorization.

diff --git a/dev/limits/index.html b/dev/limits/index.html index d1518603..5b243f28 100644 --- a/dev/limits/index.html +++ b/dev/limits/index.html @@ -37,18 +37,18 @@ ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ) ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── - 2 0 3.82 3.16 0.62 4.03e+06 -1.09e+07 0 0.775 1 1 - 4 0 3.95 6.18 0.841 1.59e+05 -1.27e+06 0 0.673 1 1 + 2 0 3.82 3.16 0.636 4.03e+06 -1.09e+07 0 0.775 1 1 + 4 0 3.95 6.18 0.843 1.59e+05 -1.27e+06 0 0.673 1 1 8 0 4.39 6.14 1.65 2.75e+05 -1.25e+05 0 0.66 1 1 - 16 0 5.36 5.75 3.27 4.66e+05 -2.33e+03 0 0.642 1 1 - 32 0 5.93 6.66 6.55 6.97e+05 -444 0 0.594 1 1 - 64 1 6.25 3.27 13.1 4.56e+07 -16.4 0 0.693 1 1 - 128 7 7.04 3.89 26.6 7.95e+07 -17.3 0 0.648 1 1 - 256 16 7.5 3.87 52.7 1.57e+08 -16.8 6.06e-305 0.633 1 1 - 512 32 7.49 4.22 105 3.12e+08 -17.2 2.2e-25 0.622 1 1 - 1.02e+03 104 8.04 2.96 211 6.23e+08 -17.1 0.00196 0.645 1 1 - 2.05e+03 189 8.42 3.34 425 1.25e+09 -17.1 0.00111 0.621 1 1 - 4.1e+03 400 8.86 3.36 851 2.5e+09 -17.2 0.0046 0.606 1 1 + 16 0 5.36 5.75 3.28 4.66e+05 -2.33e+03 0 0.642 1 1 + 32 0 5.93 6.66 6.53 6.97e+05 -444 0 0.594 1 1 + 64 1 6.25 3.27 13.3 4.56e+07 -16.4 0 0.693 1 1 + 128 7 7.04 3.89 26.8 7.95e+07 -17.3 0 0.648 1 1 + 256 16 7.5 3.87 53 1.57e+08 -16.8 6.06e-305 0.633 1 1 + 512 32 7.49 4.22 106 3.12e+08 -17.2 2.2e-25 0.622 1 1 + 1.02e+03 104 8.04 2.96 212 6.23e+08 -17.1 0.00196 0.645 1 1 + 2.05e+03 189 8.42 3.34 428 1.25e+09 -17.1 0.00111 0.621 1 1 + 4.1e+03 400 8.86 3.36 857 2.5e+09 -17.2 0.0046 0.606 1 1 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Plot the marginal mass vs. semi-major axis posterior with contours using PairPlots.jl:

pairplot(
     PairPlots.Series(
         (;
@@ -231,4 +231,4 @@
         :sma=>"log₂ Semi-major axis [au]",
         :mass=>"log₂ Mass [Mⱼᵤₚ]"
     )
-)
Example block output +)Example block output diff --git a/dev/loading-saving/index.html b/dev/loading-saving/index.html index dca825fe..df7b595a 100644 --- a/dev/loading-saving/index.html +++ b/dev/loading-saving/index.html @@ -13,4 +13,4 @@ Arrow.write("chains.arrow", tbl) # or df

You can also convert a chain object to general Array which you can save in any format you wish:

arr = Array(chains)

Saving and Restoring Models

We recommend that you save each model as a script that generates the model, e.g. in a julia file called model-systemname.jl.

For convenience, it is also possible to save and restore the full model. This is not garuanteed to work across Julia versions or between computers, but is very fast for interactive work etc.

Saving model:

using Serialization
 serialize("mymodel-systemname.jls", model)

Restoring model:

using Octofitter # must load all previously used dependencies
 using Serialization 
-model = deserialize("mymodel-systemname.jls")
+model = deserialize("mymodel-systemname.jls") diff --git a/dev/mass-photometry/index.html b/dev/mass-photometry/index.html index 125e8210..457516a7 100644 --- a/dev/mass-photometry/index.html +++ b/dev/mass-photometry/index.html @@ -31,4 +31,4 @@ M ~ Normal(1.0, 0.1) plx ~ Normal(12, 0.01) age ~ Normal(15, 1) -end b

Here the K_band_contrast_interp you supply accepts the mass of the planet, age of the system, and temperature of the planet as input, and returns the flux at K band in the same units as your images.

+end b

Here the K_band_contrast_interp you supply accepts the mass of the planet, age of the system, and temperature of the planet as input, and returns the flux at K band in the same units as your images.

diff --git a/dev/parallel-sampling/index.html b/dev/parallel-sampling/index.html index 94529214..8eba4d7b 100644 --- a/dev/parallel-sampling/index.html +++ b/dev/parallel-sampling/index.html @@ -84,4 +84,4 @@ results = Chains(model, pt) -octocorner(model, results, small=true) +octocorner(model, results, small=true) diff --git a/dev/pma/index.html b/dev/pma/index.html index bf1ff8e9..294ebd16 100644 --- a/dev/pma/index.html +++ b/dev/pma/index.html @@ -58,19 +58,19 @@ ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ) ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── - 2 0 4.53 2.1 0.0675 1.13e+05 -2.55e+03 0 0.786 1 1 + 2 0 4.53 2.1 0.0676 1.13e+05 -2.55e+03 0 0.786 1 1 4 0 4.51 6.93 0.119 1.5e+05 -20.9 0 0.631 1 1 - 8 0 6.07 7.48 0.223 2.53e+05 -15.5 2.52e-75 0.563 1 1 - 16 0 6.51 8.29 0.433 3.76e+05 -17.6 1.53e-45 0.523 1 1 - 32 0 7.52 9.24 0.853 6.21e+05 -20.1 9.87e-08 0.459 1 1 - 64 0 7.68 6.44 1.67 5.23e+07 -24.6 0.0685 0.544 1 1 + 8 0 6.07 7.48 0.224 2.53e+05 -15.5 2.52e-75 0.563 1 1 + 16 0 6.51 8.29 0.446 3.76e+05 -17.6 1.53e-45 0.523 1 1 + 32 0 7.52 9.24 0.858 6.21e+05 -20.1 9.87e-08 0.459 1 1 + 64 0 7.68 6.44 1.7 5.23e+07 -24.6 0.0685 0.544 1 1 128 2 8.63 6.65 3.38 1.06e+08 -24.6 0.169 0.507 1 1 - 256 4 8.83 6.45 6.76 2.13e+08 -22.6 0.35 0.507 1 1 - 512 20 8.67 6.49 13.4 4.23e+08 -23.2 0.327 0.511 1 1 - 1.02e+03 51 8.78 6.23 27 8.55e+08 -22.5 0.37 0.516 1 1 - 2.05e+03 99 8.66 6.01 54.1 1.7e+09 -23 0.398 0.527 1 1 - 4.1e+03 197 8.74 6.18 108 3.42e+09 -22.5 0.394 0.519 1 1 - 8.19e+03 415 8.71 6.19 216 6.82e+09 -22.5 0.409 0.519 1 1 + 256 4 8.83 6.45 6.87 2.13e+08 -22.6 0.35 0.507 1 1 + 512 20 8.67 6.49 13.6 4.23e+08 -23.2 0.327 0.511 1 1 + 1.02e+03 51 8.78 6.23 27.3 8.55e+08 -22.5 0.37 0.516 1 1 + 2.05e+03 99 8.66 6.01 54.8 1.7e+09 -23 0.398 0.527 1 1 + 4.1e+03 197 8.74 6.18 109 3.42e+09 -22.5 0.394 0.519 1 1 + 8.19e+03 415 8.71 6.19 219 6.82e+09 -22.5 0.409 0.519 1 1 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Note that octofit_pigeons took somewhat longer to run than octofit typically does; however, as we will see, it sampled successfully from severally completely disconnected modes in the posterior. That makes it a good fit for sampling from proper motion anomaly and relative astrometry with limited orbital coverage.

Analysis

The first step is to look at the table output above generated by MCMCChains.jl. The rhat column gives a convergence measure. Each parameter should have an rhat very close to 1.000. If not, you may need to run the model for more iterations or tweak the parameterization of the model to improve sampling. The ess column gives an estimate of the effective sample size. The mean and std columns give the mean and standard deviation of each parameter.

The second table summarizes the 2.5, 25, 50, 75, and 97.5 percentiles of each parameter in the model.

Pair Plot

If we wish to examine the covariance between parameters in more detail, we can construct a pair-plot (aka. corner plot).

# Create a corner plot / pair plot.
 # We can access any property from the chain specified in Variables
 using CairoMakie: Makie
@@ -90,4 +90,4 @@
             ticks=2 .^ (0:1:6)
         )
     )
-)
Example block output +)Example block output diff --git a/dev/post-pred/index.html b/dev/post-pred/index.html index 5e156487..31d1fbc1 100644 --- a/dev/post-pred/index.html +++ b/dev/post-pred/index.html @@ -32,8 +32,8 @@ Iterations = 1:1:1000 Number of chains = 1 Samples per chain = 1000 -Wall duration = 4.96 seconds -Compute duration = 4.96 seconds +Wall duration = 4.97 seconds +Compute duration = 4.97 seconds parameters = M, plx, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_θy, b_θx, b_ω, b_Ω, b_θ, b_tp internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -128,4 +128,4 @@ Makie.scatter!(ax, [0],[0], marker='⋆', color=:black, markersize=20,label="") Makie.xlims!(400,-700) Makie.ylims!(-200,200) -figExample block output

Looks like a great match to the data! Notice how the uncertainty around the middle point is lower than the ends. That's because the orbit's posterior location at that epoch is also constrained by the surrounding data points. We can know the location of the planet in hindsight better than we could measure it!

You can follow this same procedure for any kind of data modelled with Octofitter.

+figExample block output

Looks like a great match to the data! Notice how the uncertainty around the middle point is lower than the ends. That's because the orbit's posterior location at that epoch is also constrained by the surrounding data points. We can know the location of the planet in hindsight better than we could measure it!

You can follow this same procedure for any kind of data modelled with Octofitter.

diff --git a/dev/prior-pred/index.html b/dev/prior-pred/index.html index 0b172fc5..7e09bdc4 100644 --- a/dev/prior-pred/index.html +++ b/dev/prior-pred/index.html @@ -56,4 +56,4 @@ ) end -figExample block output

The heavy black crosses are our actual data, while the colored ones are simulations drawn from our priors. Notice that our real data lies at a greater separation than most draws from the prior? That might mean the priors could be tweaked.

+figExample block output

The heavy black crosses are our actual data, while the colored ones are simulations drawn from our priors. Notice that our real data lies at a greater separation than most draws from the prior? That might mean the priors could be tweaked.

diff --git a/dev/priors/index.html b/dev/priors/index.html index 6cc020d6..420164b6 100644 --- a/dev/priors/index.html +++ b/dev/priors/index.html @@ -22,8 +22,8 @@ Iterations = 1:1:1000 Number of chains = 1 Samples per chain = 1000 -Wall duration = 1.04 seconds -Compute duration = 1.04 seconds +Wall duration = 1.06 seconds +Compute duration = 1.06 seconds parameters = M, plx, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_θy, b_θx, b_ω, b_Ω, b_θ, b_tp internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -93,4 +93,4 @@ @system System1 begin plx ~ Normal(21.219, 0.060) M ~ truncated(Normal(1.1, 0.2),lower=0.1) -end b +end b diff --git a/dev/python/index.html b/dev/python/index.html index f7db846e..9ef1870e 100644 --- a/dev/python/index.html +++ b/dev/python/index.html @@ -1,2 +1,2 @@ -Using Python · Octofitter.jl

Calling from Python

This page provides some guidance on how Octofitter can be used from Python. Our general recomendation is to download Julia and copy-paste the examples as needed, but there may be cases where it useful to embed Octofitter within a larger Python project or pipeline.

In those cases, you might consider using octofitterpy. This python package uses juliacall.py to make some Octofitter functionality available in python.

Besides the model definition, most functions can be used the same in Python as in Julia. This notebook provides some examples translated into Python.

See the octofitterpy site for installation instructions.

+Using Python · Octofitter.jl

Calling from Python

This page provides some guidance on how Octofitter can be used from Python. Our general recomendation is to download Julia and copy-paste the examples as needed, but there may be cases where it useful to embed Octofitter within a larger Python project or pipeline.

In those cases, you might consider using octofitterpy. This python package uses juliacall.py to make some Octofitter functionality available in python.

Besides the model definition, most functions can be used the same in Python as in Julia. This notebook provides some examples translated into Python.

See the octofitterpy site for installation instructions.

diff --git a/dev/quick-start/index.html b/dev/quick-start/index.html index fb875c18..1c262a06 100644 --- a/dev/quick-start/index.html +++ b/dev/quick-start/index.html @@ -1,12 +1,12 @@ -Quick Start · Octofitter.jl

Quick Start (@id quick-start)

This guide introduces the key concepts in Octofitter:

  • Likelihood objects to hold your data
  • @planet and @system models to specify variables, priors, and system architecture
  • Sampling from the posterior using MCMC
  • Plotting the results
  • Saving the chain

For installation instructions, see Installation.

Example: Fit a Single Planet Orbit

Load the required packages:

using Octofitter, Distributions, CairoMakie, PairPlots

Create a PlanetRelAstromLikelihood object containing your observational data. In this case its the position of the planet relative to the star, but many other kinds of data are supported:

astrom = PlanetRelAstromLikelihood(
+Quick Start · Octofitter.jl

Quick Start (@id quick-start)

This guide introduces the key concepts in Octofitter:

  • Likelihood objects to hold your data
  • @planet and @system models to specify variables, priors, and system architecture
  • Sampling from the posterior using MCMC
  • Plotting the results
  • Saving the chain

For installation instructions, see Installation.

Example: Fit a Single Planet Orbit

Load the required packages:

using Octofitter, Distributions, CairoMakie, PairPlots

Create a PlanetRelAstromLikelihood object containing your observational data. In this case its the position of the planet relative to the star, but many other kinds of data are supported:

astrom = PlanetRelAstromLikelihood(Table(
     epoch = [50000, 50120, 50240],      # Dates in MJD
     ra = [-505.7, -502.5, -498.2],      # [mas] East positive
     dec = [-66.9, -37.4, -7.9],         # [mas] North positive
     σ_ra = [10.0, 10.0, 10.0],          # [mas] Uncertainties
     σ_dec = [10.0, 10.0, 10.0],         # [mas] Uncertainties
     cor = [0.0, 0.0, 0.0]               # RA/Dec correlations
-)

Define a planet model with orbital elements and their prior distributions:

@planet b Visual{KepOrbit} begin
+))

Define a planet model with orbital elements and their prior distributions:

@planet b Visual{KepOrbit} begin
     a ~ Uniform(0, 100)        # Semi-major axis [AU]
     e ~ Uniform(0.0, 0.5)      # Eccentricity  
     i ~ Sine()                 # Inclination [rad]
@@ -22,4 +22,4 @@
 end b

That there are many different orbit parameterizations, each requiring different of parameters names. The KepOrbit is a full 3D keplerian orbit with Campbell parameters. Visual means that we have a defined parallax distance plx that can map separations in AU to arcseconds.

Compile the model into efficient sampling code:

model = Octofitter.LogDensityModel(HD1234)

Sample from the posterior using Hamiltonian Monte Carlo (see Samplers for other options):

chain = octofit(model)

Visualize the results with orbit plots and a corner plot:

octoplot(model, chain)     # Plot orbits and data
 octocorner(model, chain)   # Corner plot of posterior

Save the results to a FITS file (see Loading and Saving Data for other formats):

savechain("output.fits", chain)

Working with Dates

Convert dates to and from Modified Julian Days using these helper functions:

mjd("2020-01-01")     # Date string to MJD
 years2mjd(2020.0)     # Decimal year to MJD
-mjd2date(50000)       # MJD to date

Next Steps

See the Tutorials section for complete examples.

+mjd2date(50000) # MJD to date

Next Steps

See the Tutorials section for complete examples.

diff --git a/dev/rel-astrom-obs/index.html b/dev/rel-astrom-obs/index.html index 40556eff..a867c3a0 100644 --- a/dev/rel-astrom-obs/index.html +++ b/dev/rel-astrom-obs/index.html @@ -87,4 +87,4 @@ b_a 10.5284 13.5713 16.8268 20.4353 29.5862 b_θ -1.5102 -1.5034 -1.5001 -1.4963 -1.4895 b_tp 872.8307 27370.7056 48895.2526 50191.1816 50401.4851 -
octoplot(model, results_obspri)
Example block output

Compare this with the previous fit using uniform priors:

octoplot(model, results_unif_pri)
Example block output

We can compare the results in a corner plot:

octocorner(model,results_unif_pri,results_obspri,small=true)
Example block output +
octoplot(model, results_obspri)
Example block output

Compare this with the previous fit using uniform priors:

octoplot(model, results_unif_pri)
Example block output

We can compare the results in a corner plot:

octocorner(model,results_unif_pri,results_obspri,small=true)
Example block output diff --git a/dev/rel-astrom/index.html b/dev/rel-astrom/index.html index 3d877819..08640a87 100644 --- a/dev/rel-astrom/index.html +++ b/dev/rel-astrom/index.html @@ -176,4 +176,4 @@

As an additional convergence test.

Analysis

As a first pass, let's plot a sample of orbits drawn from the posterior. The function octoplot is a conveninient way to generate a 9-panel plot of velocities and position:

using CairoMakie
 octoplot(model,chain)
Example block output

This function draws orbits from the posterior and displays them in a plot. Any astrometry points are overplotted.

Pair Plot

A very useful visualization of our results is a pair-plot, or corner plot. We can use the octocorner function and our PairPlots.jl package for this purpose:

using CairoMakie
 using PairPlots
-octocorner(model, chain, small=true)
Example block output

Remove small=true to display all variables.

In this case, the sampler was able to resolve the complicated degeneracies between eccentricity, the longitude of the ascending node, and argument of periapsis.

Saving your chain

Variables can be retrieved from the chains using the following sytnax: sma_planet_b = chain["b_a",:,:]. The first index is a string or symbol giving the name of the variable in the model. Planet variables are prepended by the name of the planet and an underscore.

You can save your chain in FITS table format by running:

Octofitter.savechain("mychain.fits", chain)

You can load it back via:

chain = Octofitter.loadchain("mychain.fits")
+octocorner(model, chain, small=true)Example block output

Remove small=true to display all variables.

In this case, the sampler was able to resolve the complicated degeneracies between eccentricity, the longitude of the ascending node, and argument of periapsis.

Saving your chain

Variables can be retrieved from the chains using the following sytnax: sma_planet_b = chain["b_a",:,:]. The first index is a string or symbol giving the name of the variable in the model. Planet variables are prepended by the name of the planet and an underscore.

You can save your chain in FITS table format by running:

Octofitter.savechain("mychain.fits", chain)

You can load it back via:

chain = Octofitter.loadchain("mychain.fits")
diff --git a/dev/rv-1/index.html b/dev/rv-1/index.html index cdfcf5fc..2bbef0a3 100644 --- a/dev/rv-1/index.html +++ b/dev/rv-1/index.html @@ -155,8 +155,8 @@ Iterations = 1:1:1000 Number of chains = 1 Samples per chain = 1000 -Wall duration = 29.81 seconds -Compute duration = 29.81 seconds +Wall duration = 30.42 seconds +Compute duration = 30.42 seconds parameters = M, rv0_harps, rv0_pfs, jitter_harps, jitter_pfs, b_P, b_τy, b_τx, b_mass, b_e, b_ω, b_a, b_τ, b_tp internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -203,4 +203,4 @@ fig = Octofitter.rvpostplot(model, chain) # saved to "k2_132-rvpostplot.png"Example block output

We can also plot many number of samples from the posterior:

using CairoMakie: Makie
 octoplot(model, chain)
Example block output

And create a corner plot:

using PairPlots
 using CairoMakie: Makie
-octocorner(model, chain)
Example block output

This example continues in Fit Gaussian Process.

+octocorner(model, chain)Example block output

This example continues in Fit Gaussian Process.

diff --git a/dev/rv-gp/index.html b/dev/rv-gp/index.html index 6687fc55..1a378901 100644 --- a/dev/rv-gp/index.html +++ b/dev/rv-gp/index.html @@ -164,8 +164,8 @@ Iterations = 1:1:100 Number of chains = 1 Samples per chain = 100 -Wall duration = 373.55 seconds -Compute duration = 373.55 seconds +Wall duration = 377.54 seconds +Compute duration = 377.54 seconds parameters = M, rv0_harps, rv0_pfs, jitter_harps, jitter_pfs, gp_η₁, gp_η₂, gp_η₃, gp_η₄, b_P, b_mass, b_τy, b_τx, b_e, b_ω, b_a, b_τ, b_tp internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error @@ -224,4 +224,4 @@ figscale=1.5, # make it larger alpha=0.05, # make each sample more transparent colormap="#0072b2", -) # saved to "k2_132-plot-grid.png"Example block output

Corner plot:

octocorner(model, chain, small=true) # saved to "k2_132-pairplot-small.png"
Example block output +) # saved to "k2_132-plot-grid.png"Example block output

Corner plot:

octocorner(model, chain, small=true) # saved to "k2_132-pairplot-small.png"
Example block output diff --git a/dev/rv-multi-planet/index.html b/dev/rv-multi-planet/index.html index acf9077a..86b3bfe6 100644 --- a/dev/rv-multi-planet/index.html +++ b/dev/rv-multi-planet/index.html @@ -121,4 +121,4 @@ Z2 = pt_2p.shared.reports.summary.stepping_stone[end] Z3 = pt_2p_v2.shared.reports.summary.stepping_stone[end] -Z1, Z2, Z3
(-917.2283402516634, -796.4953000419328, -793.760128696546)

As a final treat, let's animate the orbit plots. All the previous images were visualizing a single posterior draw. In this animation, we'll loop over many different samples:

Octofitter.rvpostplot_animated(model_2p_v2, results_2p_v2)
"rv-posterior.mp4"