Merge branch 'main' of http://github.com/sefffal/Octofitter.jl

docs/make.jl

@@ -5,7 +5,7 @@ pkg"dev .. ../OctofitterRadialVelocity ../OctofitterImages ../OctofitterWhereist
 Pkg.instantiate()
 Pkg.precompile()
 
-using Documenter, Octofitter
+using Documenter, Octofitter, OctofitterRadialVelocity
 
 
 makedocs(

docs/src/api.md

@@ -8,6 +8,7 @@ AstrometryLikelihood
 RadialVelocityLikelihood              
 PhotometryLikelihood
 HGCALikelihood
+ObsPriorAstromONeil2019
 Sine
 Planet
 System

ext/PlotsExt.jl

@@ -409,7 +409,7 @@ function timeplot!(
                 end
 
                 if hasproperty(obs.table, :inst_idx)
-                    idxes = length(unique(obs.table.inst_idx))
+                    idxes = maximum(obs.table.inst_idx)
                 else
                     idxes = 1
                 end

src/Octofitter.jl

@@ -42,6 +42,7 @@ include("likelihoods/system.jl")
 include("likelihoods/relative-astrometry.jl")
 include("likelihoods/photometry.jl")
 include("likelihoods/hgca.jl")
+include("likelihoods/observable.jl")
 
 include("sampling.jl")
 
@@ -72,7 +73,8 @@ function __init__()
 
         File size: 19MiB
         """,
-        "http://physics.ucsb.edu/~tbrandt/HGCA_vEDR3.fits",
+        # "http://physics.ucsb.edu/~tbrandt/HGCA_vEDR3.fits", # file is now 404
+        "https://raw.githubusercontent.com/t-brandt/orvara/master/HGCA_vEDR3.fits",
         "23684d583baaa236775108b360c650e79770a695e16914b1201f290c1826065c"
     ))
 
@@ -98,6 +100,6 @@ function __init__()
 end
 
 if VERSION != v"1.10.0-beta1"
-    include("precompile.jl")
+    # include("precompile.jl")
 end
 end

src/likelihoods/observable.jl

@@ -0,0 +1,59 @@
+#=
+This file implements the observable-based priors of
+O'Neil 2019
+"Improving Orbit Estimates for Incomplete Orbits with a New Approach 
+to Priors: with Applications from Black Holes to Planets"
+
+Currently, only the prior for relative astrometry and not radial velocity 
+is implemented.
+=#
+
+
+"""
+    ObsPriorAstromONeil2019(astrometry_table)
+
+Given a table of astrometry (in same format as `AstrometryLikelihood`,
+but only `epoch` is required), 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 orbital parameters.
+
+A Period parameter (in years) must exist in the planet model as `P`,
+in addition to the usual semi-major axis `a` (au). Typically
+this can be provided as:
+```julia
+@planet b Visual{KepOrbit} begin
+    P ~ Uniform(0.001, 1000)
+    a = cbrt(system.M * b.P^2)
+    ...
+end AstrometryLikelihood(astrom_table) ObsPriorAstromONeil2019(astrom_table)
+```
+"""
+struct ObsPriorAstromONeil2019{TTable<:Table} <: Octofitter.AbstractLikelihood
+	table::TTable
+	function ObsPriorAstromONeil2019(observations...)
+		table = Table(observations...)
+		return new{typeof(table)}(table)
+	end
+end
+export ObsPriorAstromONeil2019
+
+function Octofitter.ln_like(like::ObsPriorAstromONeil2019, θ_planet, orbit,)
+
+    jac = 0.0
+    for j in 1:length(like.table.epoch)
+        current_epoch = like.table.epoch[j]
+        M = meananom(orbit, current_epoch)
+        E = eccanom(orbit, current_epoch)
+        jac += abs(3M*(
+            eccentricity(orbit)+cos(E)
+        ) + 2*(-2+eccentricity(orbit)^2+eccentricity(orbit)*cos(E)) *sin(E))
+    end
+
+    P = period(orbit)/365.25
+    jac *= P^(1/3) / sqrt(1-eccentricity(orbit)^2);
+
+    return -2log(jac);
+end

src/likelihoods/system.jl

@@ -131,20 +131,24 @@ function make_ln_like(system::System, θ_system)
 
     planet_exprs = Expr[]
     planet_keys = Symbol[]
+    j = 0
     for i in eachindex(system.planets)
         planet = system.planets[i]
         OrbitType = _planet_orbit_type(planet)
         θ_planet = θ_system.planets[i]
         key = Symbol("planet_$i")
 
         likelihood_exprs = map(eachindex(planet.observations)) do like
-            :(
-                ll += ln_like(
+            expr = :(
+                $(Symbol("ll$(j+1)")) = $(Symbol("ll$j")) + ln_like(
                     system.planets[$(Meta.quot(i))].observations[$like],
+                    # $(system.planets[i].observations[like]),
                     θ_system.planets.$i,
                     $(key)
                 )
             )
+            j+=1
+            return expr
         end
 
         planet_contruction = quote
@@ -159,12 +163,14 @@ function make_ln_like(system::System, θ_system)
     sys_exprs = map(system.observations) do like
         # Provide the number of observations as a compile time constant
         L = Val(length(like.table))
-        :(ll += ln_like($like, θ_system, elems, $L))
+        expr = :( $(Symbol("ll$(j+1)")) = $(Symbol("ll$j")) +  ln_like($like, θ_system, elems, $L))
+        j+=1
+        return expr
     end
 
     return @RuntimeGeneratedFunction(:(function (system::System, θ_system)
 
-        ll = zero(first(θ_system))
+        ll0 = zero(first(θ_system))
 
         # Construct all orbit elements and evaluate all their individual observation likelihoods
         $(planet_exprs...)
@@ -174,7 +180,7 @@ function make_ln_like(system::System, θ_system)
 
         $(sys_exprs...)
 
-        return ll
+        return $(Symbol("ll$j"))
     end))