Plotting recipes (in order to clean up doc plots)

[mileslucas]

We have all of the information required to plot a dust law in the ExtinctionLaw structs, so we should be able to create a plotting recipe with RecipesBase to enable usage like

law = CCM89()
plot(law)

creating this plot recipe should be pretty easy and very similar to the recipes used in Polynomials.jl.

In terms of the interface, I think we need three key components

1. what law are we plotting
2. what wavelengths are we plotting
3. Are we plotting wavenumbers (inverse micron, same effective shape as frequency) or wavelengths.

Here's what I think would be easiest to accomplish that interface

1. Trivial
2. By default, use the bounds of the law. Add a recipe with bounds f(::ExtinctionLaw, a, b) and a recipe with an iterable f(::ExtincitonLaw, X)
3. Let's add a keyword argument "unit". Default is unit="invum". Other string options include unit="angstrom", unit="micron", etc. Really we should just use Unitful for parsing these, except for the "invum" case since that won't show up in unitful despite being (AFAIK) a de facto unit.

The onus for these recipes is to stop doing what we're doing with the plot docs, which is writing scripts that must be evaluated locally to generate plots that are checked out into version control and then only used for docs. A much better workflow would be to have the plots made as part of the documentation step, and ideally we can show the plotting code inside the docs for full reproducibility. After making these recipes, the docs should be able to be rewritten without the plots dir, and instead with simple plot(law) calls using Documenter's @example blocks.

[icweaver]

Hey all, I've started learning about plot recipes to try and have a go at this, but I think that I am running into a problem with the bounds. It seems that evaluating at the lower wavelength limit throws an error, while evaluating at the upper limit seems fine:

> law = GCC09()
GCC09
  Rv: Float64 3.1

> bound_lower, bound_upper = bounds(law)
(909.09, 3030.3)

> law(bound_lower)
ERROR: out of bounds of GCC09, support is over (909.09, 3030.3) angstrom

> law(bound_upper)
1.8451270235624522

Is it intended for the bounds to not be inclusive on the lower end? If so, should there be some sort of tolerance we should add to the plot recipe to avoid going out of bounds? For now, I've just been doing something like this as a toy example:

@recipe function f(t::T) where T <: ExtinctionLaw
    @inline aa_to_invum(wave::Real) = 10000 / wave
    lower, upper = bounds(t)
    w = range(lower + 0.1, upper, length=1000)
    w_invum = aa_to_invum.(w)
    w_invum, t.(w)
end

[mileslucas]

Is it intended for the bounds to not be inclusive on the lower end?

No, the intention is that law.(bounds(law)...) should not error and should be non-zero (unless the law returns 0 there)

Why is this happening?

It seems like there slight inconsistencies in the bounds and the if-else branch inside GCC09

The lower bound is 909.09 Angstrom, which is 11.000011000011 using Float64. The if-else branch for GCC09 will throw an error for anything above 11 invum.

what is the solution?

Fix the inconsistency! Either we need to change the output of bounds or we need to change the upper invum bounds.

If we change the bounds, they should become

(909.0909090909091, 3030.3030303030305)

If we change the branching, it should become

(3.3000033000032998, 11.000011000011)

How do we stop this in the future?

This should be easy to wrap into a test that can be used for all the laws, I can open a PR for this briefly. Edit: oops, this for-loop already exists, it just doesn't test if the law doesn't error.

The recipe

What you've got looks pretty good! A couple notes

• You can use the aa_to_invum within the library instead of redefining
• once the bounds are fixed you won't need to alter the start point
• You'll want to play around with the default plot attributes. The plots docs explain these decently, but the gist is: always use full names (no aliases), use --> to allow overriding, use := to disallow overriding. I would start with label --> "law_name", and xlabel, ylabel.

[icweaver]

Awesome, thanks for the PR and the recipe tips! I'm looking forward to trying these out

[icweaver]

Ok, I am getting something like this now (c0a6a23).

Is there a way to just get the model name instead of DustExtinction.<model name>?

[mileslucas]

this looks like a great start! FYI for the multi-RV plots, you should be able to condense this

begin
    fig = plot()
    Rvs = [2.0, 3.1, 4.0, 5.0, 6.0]
    for Rv in Rvs
        law = F04(Rv)
        plot!(law)
    end
    fig
end

to this!

begin
    Rvs = [2.0 3.1 4.0 5.0 6.0]
    plot(F04.(Rvs))
end

Some improvements: I'd put the number of points into a kwarg N=1000 -> x = range(..., length=N) and start working on some of those options I listed on the original post!

The pluto notebook looks pretty cool, I enjoy the html outputs

[icweaver]

Thanks for the tips! Were you thinking something more like this (8b2506c)?

haha, glad you like the notebook. btw, the latest version on master uses JuliaMono by default now!

[icweaver]

Hi Ian from the past. Looking at this with fresh eyes, I'm inclined to follow DimensionalData.jl's lead of deprecating Plots.jl support in favor of Makie.jl and their recipes system. Do other folks have thoughts one way or the other on this?

[cgarling]

I'm a matplotlib guy and have minimal experience with Plots.jl and none with Makie so I'm not clear on what the benefits / drawbacks of the two systems are. If you're volunteering to update the plots in the docs I'm happy for you to choose.

The only opinions I would offer are

1. It would be nice if we could generate the plots on-the-fly during the documentation build rather than checking the images into the repo.
2. I'd support adding a figure showing all / most of the supported extinction laws with a fixed R_V (probably 3.1) to quickly visualize how they are different / what their relevant wavelength ranges are.

[icweaver]

bet.

Sounds like you and Miles think a lot alike, haha.

I added a quick notebook here to start noodling around with some potential designs. Here's a snippet of the current impl. I've got standing:

function Makie.convert_arguments(P::PointBased, law::ExtinctionLaw) 
    aa = range(bounds(law)...; length=1_000)
    m = map(law, aa)
    invum = map(aa_to_invum, aa)
    return convert_arguments(P, invum, m) # (Plot type, x, y)
end

lplot(law::Union{CCM89, OD94, CAL00, GCC09, VCG04, FM90}; args...) = lines(
    law;
    axis = (;
        xlabel = rich("x [μm", superscript("-1"), "]"),
        ylabel = "E(B-V)",
    ),
    args...,
)

<and friends for composing other plots together>

Makie.jl primitives working

and composite plots too

Comments welcome! Can hopefully start getting unit support and more customization going after an initial design is firmed up

[cgarling]

Looks like a good start!

Is there a reason you are dispatching on long unions (e.g., lplot(law::Union{CCM89, OD94...}) rather than something more general like lplot(law::Type{<:DustExtinction.ExtinctionLaw}; args...)? It looks like you are labelling the y-axes differently so does this have to do with the units/definitions of the laws?

I'm used to seeing reddening laws expressed as A(x) / A_V as you have for F04, F19, and M14. Are the other laws actually defined as something different (maybe E(x-V)/(E(B-V)))? This confuses me (see also my issue #43) as this does not necessarily seem to be the case -- for example, our CCM89 is documented as returning "E(B-V) in magnitudes at the given wavelength relative to the extinction at 5494.5 Å" which sounds to me like E(x-V)/(E(B-V)) but in fact returns the same values as the CCM89 law implemented in the python dust_extinction package (see here for their plots) which documents their CCM89 as returning A(x) / A_V. This is related to my issue #43 -- I think it would be good to get this cleared up so we know what our y-axis labels on these plots should be.

[icweaver]

You read my mind. Yep, those unwieldy unions were just my clumsy way of trying to follow the current spec of the y-axis label units that we have

Dust stuff is far outside of my field, so I did't want to assume, but I'm all for unifying things if that seems reasonable to you!

[icweaver]

Mkay, the rest of the plots should be up now. I settled on the following spellings:

lplot => generic "law" plot
mplot => mixture model plots
dplot => dust map plot

Not super tied to these names, and it might even be better to handle this in the type system with a MixtureModel <: ExtinctionLaw and ColorModel <: ExtinctionLaw or something, but idk

[cgarling]

I think this looks right but it's not like I've gone through every model to check whether they are defined as A(x) / A_V or E(x-V)/(E(B-V)). Doing a careful check and updating docstrings should fall under #43.

This is personal preference but I would like to see the text labels a bit larger and fully black for better legibility -- they are a bit difficult to read for me ATM. Other than that I think they look good 👍

[icweaver]

Sounds good

Thanks for that feedback! Do you have a preference between one of the predefined themes here? We could also just roll our own. Personally, I like plots that don't try and mix different fonts together

[cgarling]

The default theme looks fine on that page -- did you use theme_light on your notebook? I find theme_light harder to read

[icweaver]

Sounds good. Yep, will switch it over to default. I'm also super rusty of doc stuff, so I will make sure to take some time to re-familiarize myself with best practices used here before opening a PR 👍🏾

[cgarling]

I've only done it once but my workflow for building plots and including them in docs was

1. Make a standalone source file that can be included at the top of docs/make.jl (before calling makedocs) that generates the figures and saves them somewhere
2. The present working directory (pwd) is changed to docs/build during Documenter.makedocs so paths referencing figures must be set relative to that. I usually find it easiest to move or copy the generated plots from their original locations into the build folder which you can accomplish with a hidden @example block in a documentation markdown .md file, e.g. here.
3. Once the figures are in the build folder, you can reference them like this.

Moving the figures from where they were generated into the build folder is nice so it doesn't leave plot figures scattered in the source directory when you build the docs locally rather than CI. Also docs/build is excluded from git tracking in the .gitignore folder so it also prevents the figures from being checked into CI.

There might be an easier way to manage the file paths (i.e., the difference between where the plots are generated and where they are included) but I don't know what it is. I wonder if you can just do lplot(law) in a hidden @example block and have it appear inline in the documentation?

[icweaver]

Alright, have a first crack at this up #58. It's a fairly substantial change, but hopefully for the better. Will continue discussion about it over there