r-quantities / units

Measurement units for R
https://r-quantities.github.io/units
175 stars 28 forks source link

Arithmetic ops with offset units (e.g. temperature) should issue a warning #253

Open jkeuskamp opened 3 years ago

jkeuskamp commented 3 years ago

I have encountered an issue when multiplying one temperature unit with another:

library(units)
a= set_units(20, celsius)
b=set_units(a, fahrenheit)
c=set_units(a, K)
d=set_units(2, J/K)
e=set_units(1, J/K)

I would expect that multiplying a,b, or c with d would render equal results.

Instead I find:

d a = 296.9278 [J] 313.15 [J] d b = 552.3 [J] 330.9278 [J] d * c = 296.9278 [J] 586.3 [J]

interestingly, results are correct when setting the relation to 1:

e=set_units(1, J/K) e a = 276.15 [J] 313.15 [J] e b = 276.15 [J] 313.15 [J] e * c = 276.15 [J] 313.15 [J]

Enchufa2 commented 3 years ago

Version of the package? I obtain different results.

jkeuskamp commented 3 years ago

0.6-7 installed from CRAN

jkeuskamp commented 3 years ago

@Enchufa2 in a simpler example I get the following curious result:

a = set_units(20,celsius); b = set_units(a, K)
a/b == 273.2182 [1] 
b/a == -258.4925

Do you get correct results ?

Enchufa2 commented 3 years ago

This is what you should obtain with current version, and these are correct results:

library(units)
#> udunits system database from /usr/share/udunits

a = set_units(20, celsius)
b = set_units(a, fahrenheit)
c = set_units(a, K)
d = set_units(2, J/K)

d * a
#> 313.15 [J]
d * b
#> 330.9278 [J]
d * c
#> 586.3 [J]

Why do I say these are correct? Well, you must have in mind that Celsius and Fahrenheit are relative units, which basically means that arithmetic operations with them are physically meaningless. For example, how much is 0 °C + 0 °C? As soon as you operate with Celsius or Fahrenheit, you are implicitly assuming that these are temperature differences, not temperatures. So a difference of 0 °C + a difference of 0 °C equals a difference of 0 °C, as units will tell you. This is arithmetically correct and makes physical sense. Above, you are operating with absolute (K) and implicitly relative temperatures (Celsius, Fahrenheit respectively) at the same time, which produces arithmetically correct results, but without any physical meaning.

The units package can ensure that results 1) are arithmetically correct, and 2) satisfy dimensional analysis, but could never ensure that results are physically meaningful. That's the user's job: if you want to operate with absolute temperatures (in other words, with energy per degree of freedom), then you have to convert everything to Kelvin first.

jkeuskamp commented 3 years ago

@Enchufa2 While I see your point, the result is highly counter-intuitive to me and does not match the behaviour described in the description of the package. In the abstract it says:

'When used in expression, it automatically converts units, and simplifies units of results when possible; in case of incompatible units, errors are raised.'

Celsius and Fahrenheit are temperature differences with respect to a different datum and scale. I would therefore expect that the units are not cancelled out upon division, even though the result is indeed dimensionless.

When evaluating a = set_units(20,celsius); b = set_units(a, Fahrenheit) the answer is 32.52941 [1] whereas 32.52941 [C/Fahrenheit] would more accurately describe what has just happened. Alternatively a warning could be issued that the user is trying to do something which they (probably) shouldn't be doing.

With regards to the numbers: my bad, I had posted results rendered with a different value of the initial number -- I will correct this to clarify my post

Enchufa2 commented 3 years ago

@Enchufa2 While I see your point, the result is highly counter-intuitive to me and does not match the behaviour described in the description of the package. In the abstract it says:

'When used in expression, it automatically converts units, and simplifies units of results when possible; in case of incompatible units, errors are raised.'

Conversion and simplification is possible; these units are compatible. Compatibility is defined in terms of dimensional analysis, not physical meaning.

Celsius and Fahrenheit are temperature differences with respect to a different datum and scale. I would therefore expect that the units are not cancelled out upon division, even though the result is indeed dimensionless.

A (dimensional) relationship is defined between Kelvin, Celsius and Fahrenheit, but there is no physical meaning attached to quantities. The units system has no way of knowing whether Kelvins are absolute and the others are relative to it, or Celsius ir absolute and the others are relative to it, etc. Or they all are relative to another quantity not present in the system.

When evaluating a = set_units(20,celsius); b = set_units(a, Fahrenheit) the answer is 32.52941 [1] whereas 32.52941 [C/Fahrenheit] would more accurately describe what has just happened.

That's not correct. This is a conversion, and conversion from Celsius to Fahrenheit implies a scaling factor and an offset.

Alternatively a warning could be issued that the user is trying to do something which they (probably) shouldn't be doing.

I'd say that you expect too much from the package (i.e., physical understanding), but that's not its purpose, it solely ensures arithmetical correctness between units.

Enchufa2 commented 3 years ago

When evaluating a = set_units(20,celsius); b = set_units(a, Fahrenheit) the answer is 32.52941 [1] whereas 32.52941 [C/Fahrenheit] would more accurately describe what has just happened.

Here I assume that you meant a/b, but was lost. In that case,

library(units)
#> udunits system database from /usr/share/udunits

a = set_units(20, Celsius)
b = set_units(a, fahrenheit)
units_options(simplify=FALSE)
a/b
#> 0.2941176 [°C/fahrenheit]
jkeuskamp commented 3 years ago

I see your point. Looking back I should have titled this issue 'cancelling out temperature units by division renders unexpected results'.

Nevertheless issuing a warning for division of units where x [unit_a] / set_units(x, unit_b) does not equal 1 [1] would avoid a lot of unnecessary mistakes.

Enchufa2 commented 3 years ago

Probably you meant dividing convertible units in which the conversion implies an offset (which, AFAICT right now, only happens for temperatures), otherwise, it makes little sense, because then we wouldn't be able to define a velocity, an acceleration... That would be doable. But then, why a warning for division and not for other operations? Adding up Kelvin and Celsius is IMO as bad as dividing them.

Enchufa2 commented 3 years ago

It would be doable, but looking again at the code base, it wouldn't be easy, and I don't see any efficient way. So I add it to the whishlist, but at the moment, I don't think we want to add a complex check to all arithmetic operations just to detect improper use of temperatures and issue a warning.

jkeuskamp commented 3 years ago

Thanks for your quick answers and adding issuing a warning to the wishlist. As you mention that you do not foresee this to be implemented anytime soon, I'd like to motivate a bit more why I think this issue leads to unexpected results for users such as myself:

The implicit assumption of temperature units as temperature differences is, to my knowledge, not mentioned in the package documentation. Moreover, dividing two temperature units (i.e. x [°C]/ x [K]) with this assumption renders different results in units than in e.g. WolframAlpha

in units:

a=set_units(1,celsius)
b=set_units(274.15,K)
a/b 
#> 273.1536 [1]   # <--1/274.15 + 273.15

in Wolfram Alpha:

(1 (Celsius degree difference))/(274.15 (kelvins difference)) = 0.00036476 #<-- 1/274.1529

When assuming °C to mean 'degrees Celsius of temperature' and K to mean 'kelvins of absolute temperature' Wolfram Alpha results in:

(1 (Celsius degree))/(274.1529 (kelvins))=0.00036476 [°C/K] or 
0.9963418 [-] with temperature converted to kelvins.

Note the deviation to 1 as Wolfram Alpha uses a slightly different offset for °C to K.

Enchufa2 commented 3 years ago

Interesting. So Wolfram Alpha is basically adding that layer of physical meaning, which is great. However, we depend on UDUNITS2, which doesn't have that capability, and adding something like that in the R side on top of UDUNITS2 would be too costly for all other units. For now, we could add some docs about the issue, encouraging the user to explicitly convert temperatures to the same unit before performing arithmetic operations with them.