CurrentModule = BoxLeastSquares
To install use Pkg. From the REPL, press ] to enter Pkg-mode
pkg> add BoxLeastSquaresIf you want to use the most up-to-date version of the code, check it out from main
pkg> add BoxLeastSquares#mainFirst, import the package
julia> using BoxLeastSquares
you can optionally alias the package name, too
julia> import BoxLeastSquares as BLSnow, load some data. If you don't have an estimate of the y error it will default to 1.
DocTestSetup = quote
using StableRNGs
rng = StableRNG(3351)
t = 10 .* rand(rng, 1000);
yerr = 5e-3 .* (rand(rng, 1000) .+ 1)
P = 2; t0 = 0.5; dur = 0.16; depth = 0.2;
mask = @. abs((t - t0 + 0.5P) % P - 0.5P) < 0.5dur;
y = @. ifelse(mask, 1 - depth, 1);
y .+= yerr .* randn(rng, 1000)
load_data() = t, y, yerr
end
julia> t, y, yerr = load_data(); # load data somehow
The primary interface is through the BLS method
julia> result = BLS(t, y, yerr; duration=0.16)
BLSPeriodogram
==============
input dim: 1000
output dim: 1820
period range: 0.32 - 5.014724142709022
duration range: 0.16 - 0.16
objective: likelihood
parameters
----------
index: 1633
period: 1.99930396919953
duration: 0.16
t0: 0.5001330656464655
depth: 0.19594118110109113 ± 0.0008688097746093883
snr: 225.52828804117118
log-likelihood: 27396.365214805144
to extract the parameters in a convenient named tuple use BoxLeastSquares.params
julia> BoxLeastSquares.params(result)
(index = 1633, power = 27396.365214805144, period = 1.99930396919953, duration = 0.16, t0 = 0.5001330656464655, depth = 0.19594118110109113, depth_err = 0.0008688097746093883, snr = 225.52828804117118, loglike = 27396.365214805144)
The period grid was automatically determined using autoperiod, but you can supply your own, too:
julia> periods = exp.(range(log(2) - 0.1, log(2) + 0.1, length=1000));
julia> result_fine = BLS(t, y, yerr; duration=0.12:0.01:0.20, periods=periods)
BLSPeriodogram
==============
input dim: 1000
output dim: 1000
period range: 1.809674836071919 - 2.210341836151295
duration range: 0.12 - 0.2
objective: likelihood
parameters
----------
index: 503
period: 2.001001251543549
duration: 0.168
t0: 0.4961330656464656
depth: 0.19466955969052016 ± 0.0008627202098527317
snr: 225.64622628204188
log-likelihood: 27457.6383039924
BoxLeastSquares.jl is fully compatible with Unitful.jl (although it is not a dependency of the library). For example
julia> using Unitful
julia> tu = t * u"d";
julia> results_units = BLS(tu, y, yerr; duration=(2:0.1:4)u"hr")
BLSPeriodogram
==============
input dim: 1000
output dim: 3343
period range: 0.3333333333333333 d - 4.988348864592586 d
duration range: 2.0 hr - 4.0 hr
objective: likelihood
parameters
----------
index: 2986
period: 2.0019235780121827 d
duration: 3.8000000000000003 hr
t0: 0.4916330656464656 d
depth: 0.19445716575012517 ± 0.0008692454825826517
snr: 223.70799693127577
log-likelihood: 26953.643422397385
BoxLeastSquares.BLSPeriodogram has plotting shorthands built right in- by default it will plot the period grid and the computed power
using BoxLeastSquares, Unitful, StableRNGs # hide
rng = StableRNG(3351) # hide
t = 10 .* rand(rng, 1000)u"d" # hide
yerr = 5e-3 .* (rand(rng, 1000) .+ 1) # hide
P = 2u"d"; t0 = 0.5u"d"; dur = 0.16u"d"; depth = 0.2; # hide
mask = @. abs((t - t0 + 0.5P) % P - 0.5P) < 0.5dur # hide
y = @. ifelse(mask, 1.0 - depth, 1.0) # hide
y .+= yerr .* randn(rng, 1000) # hide
results_units = BLS(t, y, yerr; duration=(2:0.1:4)u"hr") # hide
using Plots, UnitfulRecipes
plot(results_units, label="")
now let's look at how the transit model compares to the data
pars = BoxLeastSquares.params(results_units)
wrap = 0.5 * pars.period
phases = @. (mod(t - pars.t0 + wrap, pars.period) - wrap) / pars.period
inds = sortperm(phases)
model = BoxLeastSquares.model(results_units)
scatter(phases[inds], y[inds], yerr=yerr[inds],
label="data", xlabel="phase", xlim=(-0.2, 0.2), leg=:bottomright)
plot!(phases[inds], model[inds], lw=3, label="BLS model")
This code has been benchmarked against the C implementation in astropy.timeseries.bls. The C version uses OpenMP to multi-thread some parts of the core BLS algorithm, but BoxLeastSquares.jl has no threading support currently. For a fair comparison, we set OMP_NUM_THREADS to 1 for the following tests.
This first benchmark is simply the time it takes to evaluate the BLS periodogram. Periods are pre-computed using autoperiod. We simulate different sizes of data sets (x-axis) as well as different sizes of period grids (shape). This benchmark does not use units. The code can be found in bench/benchmark.jl. Here is the information for my system-
Julia Version 1.6.0
Commit f9720dc2eb* (2021-03-24 12:55 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin20.3.0)
CPU: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
OMP_NUM_THREADS = 1
JULIA_NUM_THREADS = 1using CSV, DataFrames, Plots # hide
using BoxLeastSquares # hide
benchdir = joinpath(dirname(pathof(BoxLeastSquares)), "..", "bench") # hide
results = DataFrame(CSV.File(joinpath(benchdir, "benchmark_results.csv"))) # hide
groups = groupby(results, :N_per) # hide
plot(xlabel="# data points", ylabel="time (s)") # hide
# plot main curves # hide
shapes = [:o :dtriangle :diamond] # hide
for (g, shape) in zip(groups, shapes) # hide
plot!(g.N_data, [g.t_astropy g.t_bls]; c=[1 2], shape, label="") # hide
end # hide
plot!(xscale=:log10, yscale=:log10) # hide
# create faux-legends # hide
bbox_ = bbox(0, 0, 1, 1, :bottom, :left) # hide
plot!([1 2]; c=[1 2], label=["astropy" "BoxLeastSquares.jl"], inset=(1, bbox_), # hide
bg=:transparent, border=:none, axes=false, sp=2, leg=:topleft, bgcolorlegend=:white) # hide
npers = hcat((string(k.N_per) for k in keys(groups))...) # hide
scatter!([0 0 0]; shape=shapes, c=:black, alpha=0.4, label=npers, inset=(1, bbox_), # hide
bg=:transparent, border=:none, axes=false, sp=3, ylim=(1, 2), # hide
legtitle="# periods", leg=:bottomright, legendtitlefontsize=9, bgcolorlegend=:white) # hide
If you would like to contribute, feel free to open a pull request. If you want to discuss something before contributing, head over to discussions and join or open a new topic. If you're having problems with something, open an issue.