Skip to content

Commit

Permalink
Update docs to new API
Browse files Browse the repository at this point in the history
  • Loading branch information
@sefffal
sefffal committed Apr 4, 2023
1 parent 7594c78 commit 233d69e
Show file tree
Hide file tree
Showing 8 changed files with 176 additions and 231 deletions.
99 changes: 43 additions & 56 deletions docs/src/derived.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,26 +11,22 @@ Derived variables allow you to mark certain properties as constants, reparameter
Derived variables for the system as a whole can be created as follows:

```julia
@named HD12345 = System(
Variables(
M = sys -> 1.0
plx = Normal(45., 0.02),
)
)
@system HD12345 begin
M = 1.0
plx ~ Normal(45., 0.02)
end
```
In this case, instead of including `M` as a variable in the model, we define it as a function that always returns `1.0`. This is equivalent to passing `M=1.0`.

In the following case, let's define `M` as being calculated based on another variable in the model. This is how you can do reparameterizations in Octofitter.jl
```julia
@named HD12345 = System(
Variables(
plx =Normal(45., 0.02),
logM =Normal(45., 0.02),
M = sys -> 10^sys.logM,
)
)
@system HD12345 begin
plx ~ Normal(45., 0.02)
logM ~Normal(45., 0.02)
M = 10^system.logM
end
```
We defined a new variable `logM` under priors, and then used the `Derived` block to calculate `M` from the other variables in the model.
We defined a new variable `logM` as a prior, and then calculate M from it.

In general, you can write any function you want to map from any of combination of constants and variables in the model to new variables. The only constraints are that your functions always return the same outputs for the same inputs, and are differentiable. These functions will be called in a tight loop, so you should try to make sure they are as efficient as possible.

Expand All @@ -40,18 +36,16 @@ Derived variables for an individual planet are similar, but have access to both

Here is an example of reparameterizing `e` and `a` on a planet to be logarithmic quantities:
```julia
@named b = Planet{Visual{KepOrbit}}(
Variables(
τ = Normal(0.5, 1),
ω = Normal(0.1, deg2rad(30.)),
i = Normal(0.1, deg2rad(30.)),
Ω = Normal(0.0, deg2rad(30.)),
loge = Uniform(-4, 1),
loga = Normal(1, 1)
e = (sys, pl) -> 10^pl.loge,
a = (sys, pl) -> 10^pl.loga,
),
)
@planet b Visual{KepOrbit}
τ ~ Normal(0.5, 1)
ω ~ Normal(0.1, deg2rad(30.))
i ~ Normal(0.1, deg2rad(30.))
Ω ~ Normal(0.0, deg2rad(30.))
loge ~ Uniform(-4, 1)
loga ~ Normal(1, 1)
e = 10^b.loge
a = 10^b.loga
end
```
Here `e` is defined as log-uniform, and `a` as log-normal.

Expand All @@ -61,37 +55,30 @@ Planets can have Derived variables that are calculated from variables defined on
This makes it easy to, for example, create a system of two planets that are co-planar.

```julia
@named b = Planet{Visual{KepOrbit}}(
Variables(
a = Uniform(0, 15),
e = TruncatedNormal(0, 0.1, 0, 1),
ω = Normal(0.1, deg2rad(30.)),
τ = Normal(0.5, 1),
i = (sys, pl) -> sys.i,
Ω = (sys, pl) -> sys.Ω,
),
)
@named c = Planet{Visual{KepOrbit}}(
Variables(
a = Uniform(15, 45),
e = TruncatedNormal(0, 0.1, 0, 1),
ω = Normal(0.1, deg2rad(30.)),
τ = Normal(0.5, 1),
i = (sys, pl) -> sys.i,
Ω = (sys, pl) -> sys.Ω,
),
)
@named HD12345 = System(
Variables(
plx = Normal(45., 0.02),
M = Normal(45., 0.02),
i = Normal(0.1, deg2rad(30.)),
Ω = Normal(0.0, deg2rad(30.)),
),
b, c
)
@planet b Visual{KepOrbit} begin
a ~ Uniform(0, 15)
e ~ TruncatedNormal(0, 0.1, 0, 1)
ω ~ Normal(0.1, deg2rad(30.))
τ ~ Normal(0.5, 1)
i = system.i
Ω = system.Ω
end
@planet c Visual{KepOrbit} begin
a ~ Uniform(15, 45)
e ~ TruncatedNormal(0, 0.1, 0, 1)
ω ~ Normal(0.1, deg2rad(30.))
τ ~ Normal(0.5, 1)
i = system.i
Ω = system.Ω
end
@system HD12345 begin
plx ~ Normal(45., 0.02),
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 omit those two variables from their priors, and instead just take their values from the System. This way we can enforce co-planarity between planets without e.g. rejection sampling.
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:
Expand Down
32 changes: 13 additions & 19 deletions docs/src/images.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,18 +76,15 @@ You can also provide images from multiple bands and they will be sampled indepen

Now specify the planet:
```julia
@named X = Planet{Visual{KepOrbit}}(
Variables(
a = Normal(13, 3),
e = TruncatedNormal(0.2, 0.2, 0, 1.0),
τ = Normal(0.5, 1),
ω = Normal(0.1, deg2rad(30.)),
i = Normal(0.6, deg2rad(10.)),
Ω = Normal(0.0, deg2rad(30.)),
H = Normal(3.8, 0.5)
),
image_data
)
@planet Visual{KepOrbit} begin
a ~ Normal(13, 3)
e ~ TruncatedNormal(0.2, 0.2, 0, 1.0)
τ ~ Normal(0.5, 1)
ω ~ Normal(0.1, deg2rad(30.))
i ~ Normal(0.6, deg2rad(10.))
Ω ~ Normal(0.0, deg2rad(30.))
H ~ Normal(3.8, 0.5)
end image_data
```
Note how we also provided a prior on the photometry called `H`. We can put any name we want here, as long as it's used consistently throughout the model specification.

Expand All @@ -96,13 +93,10 @@ See [Fit AstrometryLikelihood](@ref fit-astrometry) for a description of the dif

Finally, create the system and pass in the planet.
```julia
@named HD82134 = System(
Variables(
M = Normal(2.0, 0.1),
plx =Normal(45., 0.02),
),
X,
)
@system HD82134 begin
M ~ Normal(2.0, 0.1),
plx ~ Normal(45., 0.02),
end X
```

If you want to search for two or more planets in the same images, just create multiple `Planet`s and pass the same images to each. You'll need to adjust the priors in some way to prevent overlap.
Expand Down
11 changes: 9 additions & 2 deletions docs/src/loading-saving.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,10 +27,17 @@ Once loaded, you can access the underlying table using e.g. `astrom.table`.

## Saving Chains

You can save your chains using any Tables.jl compatible data package. I recommend Arrow.jl since it is a binary table format that takes up less space but is still cross compatible with other languages like Python. Another good choice is CSV.jl.
There are two ways you can save chains for later analysis. The first is a built in function that stores the chain and metadata into a FITS table. The second is converting the chain to a Table and saving it using any Tables.jl compatible package (CSV, Arrow, SQL, etc.)

### Example: Saving with metadata to FITS Table
```julia
Octofitter.savechain("mychain.fits", chain)

chain = Octofitter.loadchain(fname)
```


### Examples
### Example: Saving to CSV

Converting chain to a TypedTables.jl Table (re-exported by this package)
```julia
Expand Down
61 changes: 27 additions & 34 deletions docs/src/mass-photometry.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,46 +12,39 @@ J_band_contrast_interp(mass) = sqrt(mass)

Then, list your physical variable `mass` under priors. List the model functions under Derived. This way, the H and J band flux variables will be calculated off of your models mass parameter before getting compared to the images.
```julia
@named b = Planet{Visual{KepOrbit}}(
Variables(
a = Normal(16, 3),
e = TruncatedNormal(0.2, 0.2, 0, 1.0),
τ = Normal(0.5, 1),
ω = Normal(0.6, 0.2),
i = Normal(0.5, 0.2),
Ω = Normal(0.0, 0.2),
mass = Uniform(0, 1),
H = (sys, pl) -> H_band_contrast_interp(pl.mass),
J = (sys, pl) -> J_band_contrast_interp(pl.mass),
),
)
@planet b Visual{KepOrbit} begin
a ~ Normal(16, 3)
e ~ TruncatedNormal(0.2, 0.2, 0, 1.0)
τ ~ Normal(0.5, 1)
ω ~ Normal(0.6, 0.2)
i ~ Normal(0.5, 0.2)
Ω ~ Normal(0.0, 0.2)
mass ~ Uniform(0, 1)
H = H_band_contrast_interp(b.mass)
J = J_band_contrast_interp(b.mass)
end
```

If your model grids contain more independent variables, like age, surface gravity, etc. you can create a multi-dimensional interpolator. I recomment `ThinPlate()` from Interpolations.jl as a starting point.

This might look a little like this:
models mass parameter before getting compared to the images.
```julia
@named b = Planet{Visual{KepOrbit}}(
Variables(
a = Normal(16, 3),
e = TruncatedNormal(0.2, 0.2, 0, 1.0),
τ = Normal(0.5, 1),
ω = Normal(0.6, 0.2),
i = Normal(0.5, 0.2),
Ω = Normal(0.0, 0.2),
mass = Uniform(0, 1),
temp = Normal(1200, 500),
K = (sys, pl) -> K_band_contrast_interp(pl.mass, sys.age, pl.temp),
),
)
@named HD12345 = System(
Variables(
M = Normal(1.0, 0.1),
plx = Normal(12, 0.01),
age = Normal(15, 1),
),
b
)
@planet b Visual{KepOrbit} begin
a ~ Normal(16, 3)
e ~ TruncatedNormal(0.2, 0.2, 0, 1.0)
τ ~ Normal(0.5, 1)
ω ~ Normal(0.6, 0.2)
i ~ Normal(0.5, 0.2)
Ω ~ Normal(0.0, 0.2)
mass ~ Uniform(0, 1)
temp ~ Normal(1200, 500)
K = K_band_contrast_interp(b.mass, system.age, b.temp)
end
@system HD12345 begin
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.
32 changes: 12 additions & 20 deletions docs/src/modelling.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,18 +27,14 @@ astrom = AstrometryLikelihood(
# Or from a file:
# astrom = CSV.read("mydata.csv", AstrometryLikelihood)

@named B = Planet{Visual{KepOrbit}}(
Variables(
a = truncated(Normal(10, 4), lower=0, upper=100),
e = Uniform(0.0, 0.5),
i = Sine(),
ω = UniformCircular(),
Ω = UniformCircular(),
τ = UniformCircular(1.0),

),
astrom
)
@planet Visual{KepOrbit} begin
a ~ truncated(Normal(10, 4), lower=0, upper=100)
e ~ Uniform(0.0, 0.5)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
τ ~ UniformCircular(1.0)
end astrom
```

There's a lot going on here, so let's break it down.
Expand Down Expand Up @@ -73,13 +69,10 @@ If you have many observations you may prefer to load them from a file or databas
A system represents a host star with one or more planets. Properties of the whole system are specified here, like parallax distance and mass of the star. This is also where you will supply data like images and astrometric acceleration in later tutorials, since those don't belong to any planet in particular.

```julia
@named HD82134 = System(
Variables(
M = truncated(Normal(1.2, 0.1), lower=0),
plx = truncated(Normal(50.0, 0.02), lower=0),
),
B,
)
@system HD82134 begin
M ~ truncated(Normal(1.2, 0.1), lower=0),
plx ~ truncated(Normal(50.0, 0.02), lower=0),
end B
```

The `Variables` block works just like it does for planets. Here, the two parameters you must provide are:
Expand All @@ -89,7 +82,6 @@ The `Variables` block works just like it does for planets. Here, the two paramet
After that, just list any planets that you want orbiting the star. Here, we pass planet B.
You can name the system and planets whatever you like.

Note: the `@named` convenience macro just passes in the name as a keyword argument, e.g. `name=:HD82134`. This makes sure that the variable name matches what gets displayed in the package output, and saved a few keystrokes. (taken from ModellingToolkit.jl)

## Prepare model
We now convert our declarative model into efficient, compiled code.
Expand Down

0 comments on commit 233d69e

Please sign in to comment.