Plotting Unitful Functions

aatifsyed · October 24, 2020, 1:15am

Hi All,
Ramping up on Julia, so forgive any naivety
I’m hitting an error when trying to plot a function in a unitful way:

julia> using Unitful, UnitfulRecipes, Plots

julia> """
       The Severinghaus Oxygen Dissociation Curve
       http://www-users.med.cornell.edu/~spon/picu/calc/o2satcal.htm
       """
       severinghaus(ppO2) = (23400 * (ppO2^3 + 150 * ppO2)^-1 + 1)^-1
severinghaus

julia> """
       The Severinghaus Oxygen Dissociation Curve, but with units
       http://www-users.med.cornell.edu/~spon/picu/calc/o2satcal.htm
       """
       function severinghaus(ppO2::Unitful.Pressure)
           ppO2 = Unitful.uconvert(u"Torr", ppO2)
           sat = (23400 * (ppO2^3 + 150 * ppO2)^-1 + 1)^-1
       end
severinghaus

julia> plot(severinghaus, xlims=(1,150))
**pretty graph**

julia> plot(severinghaus, xlims=(1u"Torr",150u"Torr"))
ERROR: MethodError: no method matching adapted_grid(::Base.var"#62#63"{Base.var"#62#63"{RecipesPipeline.var"#11#12"{Symbol},typeof(severinghaus)},RecipesPipeline.var"#13#14"{Symbol}}, ::Tuple{Quantity{Float64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}},Quantity{Float64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}}})
Closest candidates are:
  adapted_grid(::Any, ::Tuple{Real,Real}; max_recursions) at /home/aatif/.julia/packages/PlotUtils/3Ttrk/src/adapted_grid.jl:14
Stacktrace:
 [1] _scaled_adapted_grid(::Function, ::Symbol, ::Symbol, ::Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}}, ::Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}}) at /home/aatif/.julia/packages/RecipesPipeline/tkFmN/src/user_recipe.jl:307
 [2] macro expansion at /home/aatif/.julia/packages/RecipesPipeline/tkFmN/src/user_recipe.jl:286 [inlined]
 [3] apply_recipe(::Dict{Symbol,Any}, ::Function, ::Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}}, ::Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}}) at /home/aatif/.julia/packages/RecipesBase/aQmWx/src/RecipesBase.jl:281
 [4] _process_userrecipes!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{typeof(severinghaus)}) at /home/aatif/.julia/packages/RecipesPipeline/tkFmN/src/user_recipe.jl:35
 [5] recipe_pipeline!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{typeof(severinghaus)}) at /home/aatif/.julia/packages/RecipesPipeline/tkFmN/src/RecipesPipeline.jl:69
 [6] _plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{typeof(severinghaus)}) at /home/aatif/.julia/packages/Plots/5K6Ue/src/plot.jl:167
 [7] plot(::Function; kw::Base.Iterators.Pairs{Symbol,Tuple{Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}},Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}}},Tuple{Symbol},NamedTuple{(:xlims,),Tuple{Tuple{Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}},Quantity{Int64,𝐌 𝐋^-1 𝐓^-2,Unitful.FreeUnits{(Torr,),𝐌 𝐋^-1 𝐓^-2,nothing}}}}}}) at /home/aatif/.julia/packages/Plots/5K6Ue/src/plot.jl:57
 [8] top-level scope at REPL[5]:1

I’ve attempted using Unitful.ustrip at various places in severinghaus, to no avail - am I trying to do something unsupported? Are there simple changes that I can make to either my code or a package that can get this working?
Thank you for any advice

Raf · October 24, 2020, 1:52am

First it would be good to simplify the problem - you should test the method first with some values before potting it!

julia> severinghaus(1u"Torr")
ERROR: DimensionError: 1 Torr^3 and 150 Torr are not dimensionally compatible

In this case the severinghaus function is kind of an empirical hack without units. Pressure cubed plus pressure, which isn’t a thing Unitful.jl can do.

But these functions exist and can be useful, so I usually do the opposite to what you have done and strip the units out in this situation. Stripping with u"Torr" makes sure its the right units fed into the equation, instead of some other Pressure type:

severinghaus(ppO2) = (23400 * (ppO2^3 + 150 * ppO2)^-1 + 1)^-1
severinghaus(ppO2::Unitful.Pressure) = severinghaus(Unitful.ustrip(u"Torr", ppO2))

plot(severinghaus, (1:1:150)u"Torr")

Also notice you need to pass a second argument to plot with unitful values - these are what is sent to the servinghouse function, plots wont actually know to convert that from just passing unitful xlims.

~~If the output actually has a unit, you can add it back at the end of the unitful aware severinghaus method~~

Edit: after discussions below, I think you should put units on all the constants instead, even though they probably aren’t in the original.

tim.holy · October 24, 2020, 7:50am

The right fix is to recognize that the constants in severinghaus must have their own units. (I’ve never understood the allergy to actually displaying the units in such conversion factors.) 150 must really be 150Torr^2, and 23400 must really be 23400/Torr^3. If you make that change, then you don’t need the uconvert, any pressure units should work just fine.

Raf · October 24, 2020, 8:07am

“Must” doesn’t always hold… sometimes units are raised to arbitrary non-integer exponents. But in this case you’re right - although probably it isn’t written with units anywhere in the literature.

tim.holy · October 24, 2020, 8:18am

sometimes units are raised to arbitrary non-integer exponents

m^(2/3) is perfectly valid (m = meters), though not if added to something with units of m^(7/2). When the units don’t work out, there is either a mistake in the calculation or an implicit assumption somewhere, e.g., “x must be measured in meters”. A common example is empirical fits (which this might be an example of), where what people do is eyeball a curve and then design a formula that matches it. What people sometimes forget is that the coefficients that emerge from this process must make sense from the standpoint of units.

An alternative to introducing units to 150 and 23400 is instead to realize that when we say “x must be measured in meters,” what we’re actually saying is that we divide by 1m before applying the formula for x. That is to say, x = xu / (1u"m"), where xu is the user-supplied value with arbitrary units of length.

Raf · October 24, 2020, 8:35am

I still don’t know what to do with the exponent on the units, say when x and z are unitless variables:

x + pressure^(1-z)

This is from an empiricle fit. In practise it’s (pressure/u"MPa")

tim.holy · October 24, 2020, 9:00am

A solution is to modify that to

x + (pressure/Torr)^(1-z)

for instance. If the first term is dimensionless, the second one must be—you simply can’t add a dimensionless number to one that has dimensions.

It’s easy to forget how many choices empirical fits involve: often the constants have units, but the constants usually come from a fitting program that can’t handle units, so people just get out their erasers and delete them. Moreover, the result also depends on scaling: for example, if you do a least squares fit of y = f(t) vs t for some measurements y taken at different values of t, you’ll come to different answers depending on how you choose to distribute your values of t. For example, whether you space them evenly on a linear or logarthmic axis is a choice that can end up having a big effect on the constants that come out of the process.

Raf · October 24, 2020, 9:11am

Of course that’s what I do But it’s essentially the same thing as using ustrip(u"Torr", pressure), and it’s never actually in the original equation in the paper. Probably because of those handy unit erasers empiricists seem to carry.

I hadn’t thought about scaling and how removing the units might hide those issues. I like how using Unitful makes these equations stand out so much. In Fortran models you just have not idea what the units are anywhere, so you wouldn’t even think about those problems.

tim.holy · October 24, 2020, 9:20am

Right, it is essentially the same thing, but it’s more user-friendly to put the units into the formula. When they’re not, it’s annoying to find something online and then have to scroll back through pages of material to find the statement “We assume the user is measuring pressure in Torr, length in meters, …”. Putting the units into the constants and formula is a simple way to make formulas self-contained.

More generally, I think it promotes good mathematics. You can often use units as a sanity check on the validity of a calculation; I can tell I’ve made a mistake if the units don’t work out. This trick is so useful that I often insert “imaginary” units into my variables and use this gut-check even when the thing I’m calculating doesn’t obviously involve physical units.