dev branch to finally complete the basic unit handling logic

[nevrome]

I created this as a development branch for @archaeothommy and me (and whoever would like to join) to work on the archchem class and its missing core features (see e.g. #12, #18, #20).

[nevrome]

What I implemented here so far:

• removing the unit names from the column names in as_archchem, after they have been parsed and turned into unit types for the respective columns.
• optionally recovering the unit names in the column names in remove_units. To write an archchem table one can now submit it to remove_units(recover_unit_names = TRUE) %>% readr::write_csv.
• assigning units to error columns: either [%] or whatever unit the column they are refering to has. To implement that I split colnames_to_constructors into parse_colnames and build_constructors. The output of parse_colnames may serve as a useful intermediate data product for diagnostic purposes.
• defining independent units for at% and wt% that get loaded with the package, using the install_unit mechanism of the units package. They are independent in the sense that they can not be converted into anything else automatically. We will have to provide conversion functions (as the one you have already implemented).

At this point I think it would be very important to test intensively. Overall this seems to be workable, but it's also a bit brittle. We need to know the edge cases, either to address them, or at least to work around them. So the three big ToDos I see right now are:

1. Write tests validating the desired behaviour of archchem tables with the code we have right now.
2. Integrate archchem with other package features e.g. atomic and oxide conversion.
3. Write a vignette that shows the features and intended uses of archchem tables. Here we can highlight what we get with this setup: A predictable structure with column-wise meta-information, automatic and safe unit conversion, unit-aware ggplot scales etc.

[archaeothommy]

Thank you so much, @nevrome, this all looks awesome. I played a bit around with it and it seems to work fine. I lacked the headspace of writing proper tests, sorry.

I would expect that the ID column is created by renaming the respective column rather than making a copy of it. Is there a reason for this behaviour?

[archaeothommy]

Found a potential bug:

I tried to connect wt% to the udunits so that except at% and oxide% all units can be handled with the package units (see commented line in zzz.R. However, this resulted in the unit of relative errors being overwritten by the unit of the column with their values. For example, when running

test_file <- system.file("extdata", "test_data_input_good.csv", package = "ASTR")
arch <- read_archchem(test_file, id_column = "Sample", context = 1:7)

results in (focus on SiO2):

<truncated>
 $ Al2O3             : Units: [wtP] num [1:14] 3.92 3.46 5.87 4.33 4.17 3.58 5.6 4.47 3.73 6.27 ...
  ..- attr(*, "archchem_class")= chr "archchem_concentration"
 $ SiO2              : Units: [wtP] num [1:14] 31.6 29.1 30.5 25.7 43.5 ...
  ..- attr(*, "archchem_class")= chr "archchem_concentration"
 $ SiO2_errSD%       : Units: [wtP] num [1:14] 4.3 2.88 5.71 2.22 4.65 3.67 3 2.87 3.44 6.12 ...
  ..- attr(*, "archchem_class")= chr "archchem_error"
 $ P2O5              : Units: [wtP] num [1:14] 0.15 NA 0.11 0.23 0.4 0.29 0.43 0.62 0.24 0.28 ...
  ..- attr(*, "archchem_class")= chr "archchem_concentration"
<truncated>

The behavior occurs here, and deactivating this if-statement in the underlying constructor function disables it (because it disables assignment of units in general).

I was not able to figure out the reason. I suspected that purrr::map2() takes the first occurrence in the look-up table, which is for the values and not the error because both have the same base. However, this also true for the desired behavior.

You have to restart the R session and reload the package to reproduce the bug. It seems that only "load_all()" does not reset the units defined on package load.

[nevrome]

I looked at your changes with 2cbbc68...a69aebe - I didn't know that there would be so much special notation to be mapped, but good that you have a solid overview here.

I'll look into the bug you describe now. I have the suspicion that I ran into this as well three weeks ago. Looks familiar.

[nevrome]

Ok - so I can reproduce your observation and I still vaguely remember that I encountered the same issue. Afaik this has something to do with overwriting the unitless unit without a custom base. That's what you effectively do with (0.01 g) / g (as g/g = 1). In this case all unitless units are one and even % in SiO2_errSD% becomes wtP.

Following this issue here, I think there must be a solution. Looking at the help file of units::install_unit there seems to be a way to define a unitless unit with a custom base:

install_unit(symbol = character(0), def = character(0), name = character(0))

def: either

• an empty definition, which defines a new base unit;
• "unitless", which defines a new dimensionless unit;
• a relationship with existing units (see details for the syntax).

So far I have left the definition of "mass_basis" empty. When I instead set it explicitly to unitless with safe_install("mass_basis", "unitless") then I seem to get the right behaviour.

<truncated>
 $ Al2O3             : Units: [wtP] num [1:14] 3.92 3.46 5.87 4.33 4.17 3.58 5.6 4.47 3.73 6.27 ...
  ..- attr(*, "archchem_class")= chr "archchem_concentration"
 $ SiO2              : Units: [wtP] num [1:14] 31.6 29.1 30.5 25.7 43.5 ...
  ..- attr(*, "archchem_class")= chr "archchem_concentration"
 $ SiO2_errSD%       : Units: [%] num [1:14] 4.3 2.88 5.71 2.22 4.65 3.67 3 2.87 3.44 6.12 ...
  ..- attr(*, "archchem_class")= chr "archchem_error"
 $ P2O5              : Units: [wtP] num [1:14] 0.15 NA 0.11 0.23 0.4 0.29 0.43 0.62 0.24 0.28 ...
  ..- attr(*, "archchem_class")= chr "archchem_concentration"
<truncated>

I'm not sure, though, if this covers what you have in mind. You write

I tried to connect wt% to the udunits so that except at% and oxide% all units can be handled with the package units

But we already handle all units with udunits, even at%. It just has this custom base which makes it sort of incompatible with other units. But that is unavoidable, right?

[nevrome]

I added my change in 3255fa2 with some more comments on this mechanism.

[archaeothommy]

But we already handle all units with udunits, even at%. It just has this custom base which makes it sort of incompatible with other units. But that is unavoidable, right?

This was imprecise phrasing from my side, sorry. I meant "udunits doing the conversion", i.e. without the need for dedicated conversion functions. And yes, they are unavoidable for these special cases.

Thanks! I will look into the details and your commit later today.

[archaeothommy]

I looked at your changes with 2cbbc68...a69aebe - I didn't know that there would be so much special notation to be mapped, but good that you have a solid overview here.

Welcome to the inconsistency of reporting... some are now outdated, some prefer writing with parentheses and others without. I tried to cover them all 🤷

[archaeothommy]

I added my change in 3255fa2 with some more comments on this mechanism.

Works as intended 😊

[archaeothommy]

The bug we observed earlier also occurs in unify_concentration_unit... I asume because the function just takes all columns that can be converted. Do you have a more elegant idea than removing columns with relative errors before unit conversion happens in the function?

[nevrome]

The bug we observed earlier also occurs in unify_concentration_unit... I asume because the function just takes all columns that can be converted. Do you have a more elegant idea than removing columns with relative errors before unit conversion happens in the function?

Well - we decided that wtP is a unitless unit like %. That means both % and wtP will be caught by the select condition in unify_concentration_unit, right?

units::ud_are_convertible(units::deparse_unit(y), "%")

That means unify_concentration_unit will adjust the units of wtP columns and % columns.

I have two questions:

1. Was it wise to define wtP as a unitless unit as I did in my last commit? It seems every conversion from and to wtP needs to be a custom defined operation in the ASTR package, not something that the units package can automatically give us. So maybe we should treat it as something entirely unconnected, like atP.
2. Is unify_concentration_unit ill-conceived? For wtP and atP conversion we decided to implement custom functions. What are the other applications of this functions that really matter? I think originally I wrote this to unify units like ppb, ppm, ppt, % etc.. I feel your understanding of what this function could and should do goes beyond that. Maybe we should delete it and leave these trivial conversions up to the user.

[archaeothommy]

2. Is unify_concentration_unit ill-conceived? For wtP and atP conversion we decided to implement custom functions. What are the other applications of this functions that really matter? I think originally I wrote this to unify units like ppb, ppm, ppt, % etc.. I feel your understanding of what this function could and should do goes beyond that. Maybe we should delete it and leave these trivial conversions up to the user.

My understanding is that wtP == %(m/m), i.e. it differs from ppm(m/m) only in the conversion factor (10^-2 instead of 10^-6). So, I think it was wise and only consequential to define is as unitless, like the others. Your function is very convenient in this regard and I think it does what you intended.

It seems every conversion from and to wtP needs to be a custom defined operation in the ASTR package, not something that the units package can automatically give us. So maybe we should treat it as something entirely unconnected, like atP.

To me, the relative errors are the odd ones out because they have a different type of total than the other values (the specific value rather than the dataset). Therefore, my suggestion is to convert relative errors to absolute errors during import and to handle by default only absolute errors. Conversion functions between relative and absolute errors should be provided, of course.
Reporting relative errors in publicatons is convenient because they are often the same across a dataset (e.g. for XRF data) but for practical purposes, you almost always convert to absolute errors anyway (e.g. for drawing error bars). Having only absolute analytical precisions also simplifies unit conversions because they can be handled the same way as the measured value when converting with e.g. unify_concentration_unit and any other conversions of SI units. It should also be the case for (the later to be implemented) error propagation in unit conversions between atP´, wtP`, because unit conversions are currently only involving multiplication by constants; but this is something that I have to look into more deeply yet.