sefffal / Octofitter.jl

Octofitter is a Julia package for performing Bayesian inference against direct images of exoplanets, relative astrometry, and astrometric acceleration of the host star.
https://sefffal.github.io/Octofitter.jl/dev
MIT License
33 stars 2 forks source link

Number of RV instrument jitter/offset parameters #26

Closed langfzac closed 3 months ago

langfzac commented 5 months ago

Is there a particular reason the number of jitter/offset terms is restricted to 10? I'm sure it's not an issue for most folks, but I have a weird use case.

sefffal commented 4 months ago

Hi @langfzac , sorry for the slow reply! It's limited to a fixed number for performance reasons but that number could be increased without any issue. Or perhaps the code could be adjusted to accept an arbitrary number without a performance impact, but that might take more work.

How many do you need? Would 25 be enough?

langfzac commented 4 months ago

Hey @sefffal, no worries, thanks for getting back! By performance do you mean you just want the offset/jitter arrays to not allocate in the likelihood function? Or is there something deeper in the code? If the former, I can whip up a PR to make the sizes arbitrary.

I'm looking for roughly 100 (again, just a weird case I want to test out).

sefffal commented 4 months ago

It is only performance, but I think there is probably a way around it. Currently we check for and load up to 10 different offsets and jitter parameters, which means the getproperty call (eg theta_system.jitter_1) is always statically known. If we did this dynamically, eg getproperty(theta_system, Symbol("jitter_$i")) it might not be type stable.

langfzac commented 4 months ago

Right, and the Symbol("jitter_$i") will allocate. How hard would it be to just collect all the jitter and offsets together, within the total parameter @NamedTuple? Like how the planet parameters are held. Then I think you could just drop-in replace the barycentric_rv_inst and jitter_inst with the proper input \theta_system fields. No new allocations and it would (should...) be type stable.

sefffal commented 4 months ago

If it's just for a one-off use case, maybe the allocation isn't a huge deal. You could give the dynamic version a try and see if it's fast enough for you.

One option to get an unlimited amount of instruments to work would be to write something that loops over theta_system. Since all the properties of the named tuple are part of the type, I think that loop could be type-stable.

A second option would be to try and define the jitter using eg. an MvNormal(Diagonal(I(100))), that is a single distribution with multiple values. I'm not sure if that works right now but could probably be made to work.

Long term I hope to refactor things so that likelihood objects can have their own variables blocks like the system and planet do now.

sefffal commented 4 months ago

Looks like the second option with multivariate distributions is tricky, actually, and would need some thought.

langfzac commented 4 months ago

Ok, this seems to work just fine for now: https://github.com/langfzac/Octofitter.jl/commit/74aedeea48d570cf553ecf53ce92c6021185edf3. It makes the type a little messy though.

sefffal commented 4 months ago

Great idea! I think that could work, but might not quite yet be type stable. I have a couple suggestions:

We could have:

offset_symbols::TOffsets
jitter_symbols::TJitters

where the struct definition is expanded like:

struct StarAbsoluteRVLikelihood{TTable<:Table,GP,TOffsets,TJitters} <: Octofitter.AbstractLikelihood

Then the function definition can lose:

    num_inst::Val{NI}=Val(length(rvlike.instrument_names))
) where {L,N,NI}

The reason is because the default value I think is actually not type stable unless already known in the system likelihood function. The default value of num_epochs is provided for convenience but is, I think, only type-stable when it's filled in by the make_ln_likelihood function.

langfzac commented 4 months ago

Yeah, that's definitely better. I can start a RP, and we can iterate on it, if this is something you want to merge into the main repo. Otherwise, I think we can close this issue and I can just use my fork. Thanks for the help!

sefffal commented 4 months ago

Great, thanks so much! A PR would definitely be appreciated.