Added F99 model #30

icweaver · 2020-07-07T06:39:53Z

Add support for Measurements
Generalize for F04

Contributing to #10

… and indexing

codecov · 2020-07-07T07:02:19Z

Codecov Report

Merging #30 into master will increase coverage by 0.09%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #30      +/-   ##
==========================================
+ Coverage   99.35%   99.45%   +0.09%     
==========================================
  Files           4        4              
  Lines         155      182      +27     
==========================================
+ Hits          154      181      +27     
  Misses          1        1

Impacted Files	Coverage Δ
src/DustExtinction.jl	`100.00% <ø> (ø)`
src/color_laws.jl	`100.00% <100.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ad510f3...2101f98. Read the comment docs.

giordano · 2020-07-07T07:32:01Z

Everything looks to be working alright except for Measurements handling. I think this might be related to the package (Dierckx.jl) I am using for interpolation?

I'm not at the computer right now, but if I remember correctly Dierckx.jl is based on an external C library, right? If so, yeah, Measurements.jl can't propagate the uncertainty through ccalls, but you can use the @uncertain macro: https://juliaphysics.github.io/Measurements.jl/stable/examples/#@uncertain-Macro-1

src/color_laws.jl

Project.toml

test/color_laws.jl

icweaver · 2020-07-07T10:18:16Z

Everything looks to be working alright except for Measurements handling. I think this might be related to the package (Dierckx.jl) I am using for interpolation?

I'm not at the computer right now, but if I remember correctly Dierckx.jl is based on an external C library, right? If so, yeah, Measurements.jl can't propagate the uncertainty through ccalls, but you can use the @uncertain macro: https://juliaphysics.github.io/Measurements.jl/stable/examples/#@uncertain-Macro-1

This is so cool! @uncertain is working like a charm for single wavelengths. For example, I can do this:

>>> @uncertain F99()(wave_unc[1])
0.08548387096774193 ± 0.0

but I think I am having trouble broadcasting this when I try:

>>> @uncertain F99().([wave_unc[1], wave_unc[2]])
ERROR: MethodError: no method matching measurement(::Tuple{Array{Measurement{Float64},1}})

for example.

giordano · 2020-07-07T10:20:46Z

but I think I am having trouble broadcasting this when I try:

See #30 (comment) 😉

icweaver · 2020-07-07T10:23:59Z

but I think I am having trouble broadcasting this when I try:

See #30 (comment) 😉

hahaha, was just about to respond. Thanks! It's like you read my mind

re: JuliaAstro#30 (comment) Co-authored-by: Mosè Giordano <giordano@users.noreply.github.com>

src/color_laws.jl

mileslucas

Thanks for the contribution- I've got some comments about the internals, but the documentation and testing looks good at first glance.

First, I think if we can make the interpolation work using Interpolations.jl instead of Dierckx.jl we'd be better off because as a pure-julia library we avoid issues like you saw with Measurements.jl as well as auto-differentiation. However, currently Interpolations only manages gridded interpolation up to Linear degree and requires an equally-spaced grid for BSplines up to Cubic. I have no issues including this now but I think we will want to note that it is not AD friendly.

Second, keep an eye out for moving calculations earlier. This can be achieved by the use of const variables as well as putting more calculations into the F99 struct. What I mean by this is anything that can be fully specified using Rv only (such as opt_axebv_y can be added as a field to the struct- e.g. F99(Rv) = F99(Rv, f99_opt_coeff(rv)) where that function calls what you've written out. This way, even something as innocuous as F99().(wavelengths) saves by doing that calculation once instead of every iterated method call.

Also, in general, using static data structures like Tuple is better for static data, although I know not all of this code follows that. It would be useful to check performance using Vector vs Tuple vs SVector from StaticArrays for the internals some day.

mileslucas · 2020-07-08T18:32:23Z

src/color_laws.jl

+    # x value above which FM90 parametrization used
+    x_cutval_uv = 10000.0 / 2700.0
+
+    # required UV points for spline interpolation
+    x_splineval_uv = 10000.0 ./ [2700.0, 2600.0]
+


if these are constants, I think we should move them outside the function definition and declare as constants like so

const x_cutval_uv = aa_to_invum(2700) const x_splineval_uv = aa_to_invum.((2700, 2600))

additionally by manually writing out a Tuple, we avoid allocating an array each method call.

mileslucas · 2020-07-08T18:36:05Z

src/color_laws.jl

+    optnir_axav_x::Vector{<:Real},
+    optnir_axav_y::Vector{<:Real},


In general we want to use abstract types for function dispatch, such as AbstractVector{<:Real} here. However, you don't really need to type any of these variables because they're constrained by their behavior already in the code- the most generic program actually has no types!

I would remove all the type specializations here, and you can callback @code_warntype _curve_F99_method(x, Rv, ...) from the REPL to see if you have type-stability, which is what matters.

mileslucas · 2020-07-08T18:37:47Z

src/color_laws.jl

+    # FM90 model and values
+    fm90_model = FM90(c1=c1, c2=c2, c3=c3, c4=c4, x0=x0, gamma=gamma)
+    # evaluate model and get results in A(x)/A(V)
+    axav_fm90 = fm90_model.(10000.0 ./ xuv) / Rv .+ 1.0 # Expects xuv in Å


style: using the @. macro makes this line simpler and avoids bugs like the missing ./ before Rv

Suggested change

axav_fm90 = fm90_model.(10000.0 ./ xuv) / Rv .+ 1.0 # Expects xuv in Å

axav_fm90 = @. fm90_model(aa_to_invum(xuv)) / Rv + 1 # Expects xuv in Å

src/color_laws.jl

mileslucas · 2020-07-08T19:11:13Z

src/color_laws.jl

+    optnir_axav_x = 10000.0 ./
+        [26500.0, 12200.0, 6000.0, 5470.0, 4670.0, 4110.0]


Similar to others- make const if possible and use tuples to avoid allocations

const optnir_axav_x = aa_to_invum.((26500, 12200, 6000, 5470, 4670, 4110))

mileslucas · 2020-07-08T19:12:22Z

src/color_laws.jl

+         0.701 + 1.0016 * Rv,
+         1.208 + 1.0032 * Rv - 0.00033 * (Rv^2),
+    ]
+    nir_axebv_y = [0.265, 0.829] * Rv / 3.1


you'll want to make use of broadcasting here-

nir_axebv_y = @. (0.264, 0.829) * Rv / 3.1

icweaver · 2020-07-08T19:57:38Z

Thanks so much for these comments, @mileslucas! I was trying to strike a balance between doing a 1-1 translation from the original astropy code and making things more Julian. Your suggestions have helped answer a lot of "best practices" questions I've had when starting with this language and I'll start implementing them ASAP

… macro, switched to Tuples

icweaver · 2020-07-09T04:53:14Z

Sorry, not sure why this latest commit failed on some platforms while working on others. The local tests I ran before committing all passed.

…99 by sheer luck

icweaver · 2020-07-09T08:56:35Z

Sorry, not sure why this latest commit failed on some platforms while working on others. The local tests I ran before committing all passed.

Ahh, didn't realize evalpoly is only for 1.4

https://docs.julialang.org/en/v1/base/math/#Base.Math.evalpoly

That's a bummer. Being able to broadcast was super useful for re-writing this

opt_axebv_y = (
    -0.426 + 1.0044 * Rv,
    -0.050 + 1.0016 * Rv,
     0.701 + 1.0016 * Rv,
     1.208 + 1.0032 * Rv - 0.00033 * (Rv^2),
)

as

opt_axebv_y = evalpoly.(Rv, opt_axebv_y_params)

and just pack all those constants into a nested Tuple outside the main function. Is there a way to broadcast with @evalpoly or should I just drop this? I thought maybe a comprehension like this

(@evalpoly Rv param... for param in opt_axebv_y_params)

might work, but I am pretty sure I am butchering the syntax here.

giordano · 2020-07-10T00:45:51Z

src/color_laws.jl

+    end
+
+    # terms depending on Rv
+    c2 = @evalpoly (1. / Rv) f99_Rv_params...


This cannot work with this macro, you need to list the elements of the vector one by one. That's why the function was introduced

Oh, just saw your comment! Yep, just realized that when I looked at Travis. Oh well =\

That's neat though that you can do exactly that in v1.4

Alrighty, everything should be switched over now. I've never had to think about things like backwards compatibility before, so this has been an extremely useful lesson for me, thanks!

icweaver · 2020-07-10T00:47:37Z

Ah, I see now it's failing because you can't splat a Tuple into @evalpoly in Julia v1.0.5

siddharthlal25 · 2020-07-10T14:47:24Z

According to the discussion here, it looks like Interpolations.jl does not have support for cubic B-splines (which is used by the python implementation), while this package does.

I was just going through Interpolations.jl and found this, is this what you were looking for?

icweaver · 2020-07-10T16:16:08Z

According to the discussion here, it looks like Interpolations.jl does not have support for cubic B-splines (which is used by the python implementation), while this package does.

I was just going through Interpolations.jl and found this, is this what you were looking for?

Right, so my understanding was that the issue was with the next part down, which I believe @mileslucas mentions here, re: cubic splines on a non-uniform grid?

I added using Interpolations.jl to my general TODO for #10. It would be awesome if we could use it in the future! Although, it looks like it is not being actively maintained right now =\

Found an open issue here!

icweaver · 2020-07-10T18:32:29Z

Following up on the spline issue, there is this package that could be useful!

Pros: It's written in pure Julia (which should make it play nice with AD?) and can do cubic splines on an irregular grid
Cons?: Those cubic splines are actually Akima splines, so our results will differ from the astropy implementation of the dust models

I have no experience using Akima splines, but it looks like it may actually perform better than traditional cubic B-splines.

I guess it comes down to: should DustExtinction.jl exactly reproduce its Python cousin's results, or should it use the "latest and greatest" technology? tbh, I can see arguments for both.

mileslucas

Things are looking pretty good, keep an eye out for re-ordering the computations like I mentioned in the previous review.

mileslucas · 2020-07-10T22:06:42Z

src/color_laws.jl

+# x value above which FM90 parametrization used
+const f99_x_cutval_uv = aa_to_invum(2700)
+
+# required UV points for spline interpolation
+const f99_x_splineval_uv = aa_to_invum.((2700, 2600))
+
+# spline points
+const f99_optnir_axav_x = aa_to_invum.((26500, 12200, 6000, 5470, 4670, 4110))
+const f99_nir_axebv_y_params = @. (0.265, 0.829) / 3.1
+
+# c1-c2 correlation terms
+const f99_c3 = 3.23
+const f99_c4 = 0.41
+const f99_x0 = 4.596
+const f99_gamma = 0.99
+


I think these would fit in better closer to the F99 definition- you'll see things are generally clustered by law already

For sure 804fb33

…src/color_laws.jl

icweaver · 2020-07-14T16:00:01Z

It looks like the CI sporadically fails when installing Julia, but otherwise I think everything should be addressed now for another review. Thanks!

mileslucas

This is looking good! Thanks again for the contribution. The only changes left are completely internal and I'm not too concerned. Unless anybody else has concerns I'll merge this later tonight and do a new release.

mileslucas · 2020-07-14T18:53:20Z

src/color_laws.jl

+        # using cubic spline anchored in UV, optical, and IR
+
+        # optical/NIR points in input x
+


not really following the comments here

Oh, whoops! That's from some leftover Python styling before I switched to thinking in terms of "kernels". Thanks for catching that, should be fixed now! (5a6b6fd)

mileslucas · 2020-07-14T19:19:52Z

Actually, if you want to go ahead and edit Project.toml and bump the version to 0.8.0 we'll be prepped for a release.

icweaver · 2020-07-14T19:31:49Z

This is looking good! Thanks again for the contribution. The only changes left are completely internal and I'm not too concerned. Unless anybody else has concerns I'll merge this later tonight and do a new release.

Awesome, thanks! I opened up a separate PR for F04 (#31) and F19 (#32) as well. I tried to keep in mind everything y'all have taught me =]

icweaver added 6 commits July 6, 2020 03:07

added F99

62ab1b0

switched to explicit conditionals and more Julian array indexing

03c6429

added figure

7cb154f

added some initial tests

498e641

added/cleaned up function signatures, removed unnecessary parentheses…

9f6433f

… and indexing

added documention and bibliography

a682f91

icweaver changed the title ~~F99~~ Adding F99 model Jul 7, 2020

icweaver changed the title ~~Adding F99 model~~ [WIP] Adding F99 model Jul 7, 2020

icweaver changed the title ~~[WIP] Adding F99 model~~ [WIP]: Adding F99 model Jul 7, 2020

icweaver mentioned this pull request Jul 7, 2020

Implement More Methods! #10

Open

10 tasks

icweaver changed the title ~~[WIP]: Adding F99 model~~ [WIP]: added F99 model Jul 7, 2020

added citation to index

490345a

siddharthlal25 reviewed Jul 7, 2020

View reviewed changes

src/color_laws.jl Show resolved Hide resolved

src/color_laws.jl Outdated Show resolved Hide resolved

combined part of model into single if-else block

1901d92

giordano reviewed Jul 7, 2020

View reviewed changes

Project.toml Show resolved Hide resolved

test/color_laws.jl Outdated Show resolved Hide resolved

icweaver and others added 3 commits July 7, 2020 06:26

switched to rand noise in test

6830206

re: JuliaAstro#30 (comment) Co-authored-by: Mosè Giordano <giordano@users.noreply.github.com>

re-ran figure to be in sync with Github side commit 6830206

3842aae

added compatibility bounds to Dierckx.jl

cbd9d40

siddharthlal25 reviewed Jul 7, 2020

View reviewed changes

src/color_laws.jl Outdated Show resolved Hide resolved

icweaver added 2 commits July 7, 2020 13:18

added more to docstring

f458a34

dropped redundant low-level bounds check and tests

d81b97f

mileslucas requested changes Jul 8, 2020

View reviewed changes

icweaver added 2 commits July 8, 2020 17:22

moved constants to top level, dropped function signatures, applied @.…

2f4eceb

… macro, switched to Tuples

put back internal bounds check

a6da64e

icweaver added 2 commits July 8, 2020 22:55

specialized constants (F04 will use similar names)

da7f9ea

moved more constants outside of function and switched to evalpoly

e08bc31

icweaver added 3 commits July 9, 2020 01:02

updated const naming conventions to be more consistent

384e9ab

Made change to algorithm shared by F04. Was working before for just F…

558893a

…99 by sheer luck

fixed typo

86e4883

reverted back to @evalpoly for Julia v1.0 compatiblity

49696a2

giordano reviewed Jul 10, 2020

View reviewed changes

dropped Tuple splatting in @evalpoly to support Julia v1

526978a

mileslucas reviewed Jul 10, 2020

View reviewed changes

icweaver added 5 commits July 10, 2020 19:55

moved consts closer to model

804fb33

generalized function shared by F99 and F04

49de26a

logically grouped more constants

209808b

dropped unneeded F99 from plotting script and changed order to match …

9bd6d14

…src/color_laws.jl

removed extra new line

9674172

icweaver requested a review from mileslucas July 14, 2020 15:58

mileslucas approved these changes Jul 14, 2020

View reviewed changes

removed leftover python style

5a6b6fd

updated DustExtinction.jl to v0.8.0

2101f98

giordano changed the title ~~[WIP]: added F99 model~~ Added F99 model Jul 14, 2020

mileslucas merged commit c30f719 into JuliaAstro:master Jul 15, 2020

icweaver deleted the F99 branch July 15, 2020 13:17

		optnir_axav_x::Vector{<:Real},
		optnir_axav_y::Vector{<:Real},

	axav_fm90 = fm90_model.(10000.0 ./ xuv) / Rv .+ 1.0 # Expects xuv in Å
	axav_fm90 = @. fm90_model(aa_to_invum(xuv)) / Rv + 1 # Expects xuv in Å

		optnir_axav_x = 10000.0 ./
		[26500.0, 12200.0, 6000.0, 5470.0, 4670.0, 4110.0]

		# using cubic spline anchored in UV, optical, and IR

		# optical/NIR points in input x

Added F99 model #30

Added F99 model #30

Conversation

icweaver commented Jul 7, 2020 • edited Loading

codecov bot commented Jul 7, 2020 • edited Loading

Codecov Report

giordano commented Jul 7, 2020

icweaver commented Jul 7, 2020

giordano commented Jul 7, 2020

icweaver commented Jul 7, 2020 • edited Loading

mileslucas left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

icweaver commented Jul 8, 2020

icweaver commented Jul 9, 2020

icweaver commented Jul 9, 2020 • edited Loading

Choose a reason for hiding this comment

icweaver Jul 10, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

icweaver commented Jul 10, 2020

siddharthlal25 commented Jul 10, 2020

icweaver commented Jul 10, 2020 • edited Loading

icweaver commented Jul 10, 2020 • edited Loading

mileslucas left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

icweaver commented Jul 14, 2020

mileslucas left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

icweaver Jul 14, 2020 • edited Loading

Choose a reason for hiding this comment

mileslucas commented Jul 14, 2020

icweaver commented Jul 14, 2020 • edited Loading

icweaver commented Jul 7, 2020 •

edited

Loading

codecov bot commented Jul 7, 2020 •

edited

Loading

icweaver commented Jul 7, 2020 •

edited

Loading

icweaver commented Jul 9, 2020 •

edited

Loading

icweaver Jul 10, 2020 •

edited

Loading

icweaver commented Jul 10, 2020 •

edited

Loading

icweaver commented Jul 10, 2020 •

edited

Loading

icweaver Jul 14, 2020 •

edited

Loading

icweaver commented Jul 14, 2020 •

edited

Loading