rounding issue in fwrite

[haohua]

Hi guys,

when i am writing some numerical data out to text file with data.table::fwrite, it lost precision, which is by design because of the 16 decimal point limit.

But, here is a confusing part I have for fwrite, two set of data, one with 0.xxxx, and the other with 5.xxxx.
The output file gives different ranking order, compared with the original data.

library(data.table)
testNum1 = c(0.87804878048780566,
             0.87804878048780477,
             0.87804878048780477)
fwrite(x = data.table(value = testNum1), 'test1.txt')
testNumRead1 = fread('test1.txt')

testNum2 = c(5.87804878048780566,
             5.87804878048780477,
             5.87804878048780477)
fwrite(x = data.table(value = testNum2), 'test2.txt')
testNumRead2 = fread('test2.txt')


frank(testNum1)
# [1] 3.0 1.5 1.5
frank(testNumRead1)
[1] 3.0 1.5 1.5

frank(testNum2)
# [1] 3.0 1.5 1.5
frank(testNumRead2)
[1] 2 2 2

can anyone help me to understand the issue here?

[MichaelChirico]

It may not look like much, but adding the 5 in front changes how fwrite handles your number because of significant figures:

data.table/src/fwrite.c

Lines 26 to 27 in da7fca8

	#define NUM_SF 15
	#define SIZE_SF 1000000000000000ULL // 10^NUM_SF

Adding a number to the left of the decimal point means the significant figure count reaches the max before the place where the decimals differ. That is, it's not a 16 decimal point limit, but a 16 significant figure limit

You could also have added a 0 in the middle of the first number for the same effect:

val = c(
  0.878048780487800566,
  0.878048780487800477,
  0.878048780487800477
)
fwrite(data.table(val1))
# val
# 0.878048780487801
# 0.878048780487801
# 0.878048780487801

The way to fix this would be to un-fix the NUM_SF parameters which are hard-coded now.

Note that scipen does not fix the issue here, because the number is truncated to significant digits before scipen is even considered.

Also note that the numbers are in fact distinct in double representation -- see https://float.exposed for the example numbers. So we could in fact write them with more precision.

leaving this as an FR for now

[MichaelChirico]

related: #2398

[niekverw]

I am running into the same issue.

Is there a way to convert the data.table such that it is consistent with what fwrite would write, then I can predict the behavior and work around it. i.e. is it possible to fwrite to memory ? or execute the rounding function from fwrite without actually using fwrite to a file?

[ben-schwen]

I am running into the same issue.

Is there a way to convert the data.table such that it is consistent with what fwrite would write, then I can predict the behavior and work around it. i.e. is it possible to fwrite to memory ? or execute the rounding function from fwrite without actually using fwrite to a file?

We already have a fast numeric to string binary function accessible via data.table:::binary.

@niekverw
You could use the code below to define the reverse operation and use to serialize/unserialize your data.

library(data.table)
library(inline)
src = "
  if (!isString(x)) error(\"x must be type 'character'\");
  int n = LENGTH(x);
  SEXP ans = PROTECT(allocVector(REALSXP, n));
  double *ansd = REAL(ans);
  for (int i=0; i<n; i++) {
    const char *str = CHAR(STRING_ELT(x, i));
    uint64_t i64 = 0;
    for (int j=0; j<68; j++) {
      if (str[j] == ' ') continue;
      i64 <<= 1;
      if (str[j] == '1') i64 |= 1;
      else if (str[j] != '0') error(\"Found invalid binary representation while converting to decimal\");
    }
    double result; // directly casting pointer breaks strict aliasing
    memcpy(&result, &i64, sizeof(double));
    ansd[i] = result;
  }
  UNPROTECT(1);
  return ans;
"
decimal = cfunction(signature(x = 'character'), body=src, includes="#include <stdint.h>")

is_bin = function(x) {
    if (!is.character(x) || nchar(x) != 68L) return(FALSE)
    tbl = table(strsplit(x, '')[[1]])
    tbl[1] == 4L && all(names(tbl) %in% c(" ", "0", "1"))
}
is.binary = function(x) all(vapply(x, is_bin, logical(1)))

dt = data.table(a = letters[1:3], b=c(0, NA_real_, 0.5), c=c(0.87804878048780566, 0.87804878048780477, 0.87804878048780477))
cols = names(which(sapply(dt, is.numeric)))
dt[, eval(cols) := lapply(.SD, data.table:::binary), .SDcols = cols]
fwrite(dt, f<-tempfile())

dt2 = fread(f)
cols2 = names(which(sapply(dt, is.binary)))
dt2[, eval(cols2) := lapply(.SD, decimal), .SDcols = cols2]
dt2
#>         a     b         c
#>    <char> <num>     <num>
#> 1:      a   0.0 0.8780488
#> 2:      b    NA 0.8780488
#> 3:      c   0.5 0.8780488

frank(dt$c)
#> [1] 3.0 1.5 1.5
frank(dt2$c)
#> [1] 3.0 1.5 1.5

[tdhock]

there are inherent limits in R real numbers (C double) which need to be considered

[dvg-p4]

data.table is currently rounding numbers more than the inherent limits of double-precision floating points though, in a way that can actively corrupt the data. I'd consider this a bug:

> data.table::fwrite(data.frame(num = 1234567890123456), "", row.names = FALSE)
num
1234567890123460

Note that base-R write.table() does not make the same mistake:

> write.table(data.frame(num = 1234567890123456), "", row.names = FALSE)
"num"
1234567890123456

[dvg-p4]

It's definitely a bit of an edge case here, since it only happens to integers that are close (on a log scale) to the maximum that doubles can precisely represent. In a programming language that actually had a built-in 64-bit integer type, using doubles for storing large precise integers would be an antipattern in the first place...Still, would be good to handle them at least as well as write.table() can.

Ah, and for the record, looks like the bit64::integer64 representation of the number is printed correctly, so if the use of bit64 is acceptable for one's use case, that should be a decent workaround.

> data.table::fwrite(data.frame(num = bit64::as.integer64(1234567890123456)), "", row.names = FALSE)
num
1234567890123456