# Plotting Unitful Functions

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:
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?

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.

1 Like

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.

6 Likes

β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.

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.

5 Likes

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")`

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.

3 Likes

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.

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.

7 Likes