Implement and use the hts_str2dbl method (by @jkbonfield)

For faster conversion of strings to floating point values.  It
speeds up the conversion of values which match the regexp
[+-]?[0-9]*[.]?[0-9]* and have at least one and no more than 15
digits in total.  Values that don't match (for example they include
an exponent) are handled by passing to strtod().

Co-Authored-By: James Bonfield <[email protected]>
Co-Authored-By: Rob Davies <[email protected]>

Author: 3 people

textutils_internal.h

@@ -1,6 +1,6 @@
 /* textutils_internal.h -- non-bioinformatics utility routines for text etc.
 
-   Copyright (C) 2016,2018,2019 Genome Research Ltd.
+   Copyright (C) 2016,2018-2020 Genome Research Ltd.
 
    Author: John Marshall <[email protected]>
 
@@ -304,6 +304,96 @@ static inline uint64_t hts_str2uint(const char *in, char **end, int bits,
     return n;
 }
 
+/// Convert a string to a double, with overflow detection
+/** @param[in]  in     Input string
+    @param[out] end    Returned end pointer
+    @param[out] failed Location of overflow flag
+    @return String value converted to a double
+
+Converts a floating point value string to a double.  The string should
+have the format [+-]?[0-9]*[.]?[0-9]* with at least one and no more than 15
+digits.  Strings that do not match (inf, nan, values with exponents) will
+be passed on to strtod() for processing.
+
+If the value is too big, the largest possible value will be returned;
+if it is too small to be represented in a double zero will be returned.
+In both cases errno will be set to ERANGE.
+
+If no characters could be converted, *failed will be set to 1.
+
+The address of the first character following the converted number will
+be stored in *end.
+
+Both end and failed must be non-NULL.
+ */
+
+static inline double hts_str2dbl(const char *in, char **end, int *failed) {
+    uint64_t n = 0;
+    int max_len = 15;
+    const unsigned char *v = (const unsigned char *) in;
+    const unsigned int ascii_zero = '0'; // Prevents conversion to signed
+    int neg = 0, point = -1;
+    double d;
+    static double D[] = {1,1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7,
+                         1e8, 1e9, 1e10,1e11,1e12,1e13,1e14,1e15,
+                         1e16,1e17,1e18,1e19,1e20};
+
+    while (isspace(*v))
+        v++;
+
+    if (*v == '-') {
+        neg = 1;
+        v++;
+    } else if (*v == '+') {
+        v++;
+    }
+
+    switch(*v) {
+    case '1': case '2': case '3': case '4':
+    case '5': case '6': case '7': case '8': case '9':
+        break;
+
+    case '0':
+        if (v[1] != 'x' && v[1] != 'X') break;
+        // else fall through (hex number)
+
+    default:
+        // Non numbers, like NaN, Inf
+        d = strtod(in, end);
+        if (*end == in)
+            *failed = 1;
+        return d;
+    }
+
+    while (*v == '0') ++v;
+
+    const unsigned char *start = v;
+
+    while (--max_len && *v>='0' && *v<='9')
+        n = n*10 + *v++ - ascii_zero;
+    if (max_len && *v == '.') {
+        point = v - start;
+        v++;
+        while (--max_len && *v>='0' && *v<='9')
+            n = n*10 + *v++ - ascii_zero;
+    }
+    if (point < 0)
+        point = v - start;
+
+    // Outside the scope of this quick and dirty parser.
+    if (!max_len || *v == 'e' || *v == 'E') {
+        d = strtod(in, end);
+        if (*end == in)
+            *failed = 1;
+        return d;
+    }
+
+    *end = (char *)v;
+    d = n / D[v - start - point];
+
+    return neg ? -d : d;
+}
+
 
 #ifdef __cplusplus
 }

vcf.c

@@ -2201,7 +2201,7 @@ static int vcf_parse_format(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v, char *p
 {
     if ( !bcf_hdr_nsamples(h) ) return 0;
 
-    static int extreme_int_warned = 0;
+    static int extreme_val_warned = 0;
     char *r, *t;
     int j, l, m, g, overflow = 0;
     khint_t k;
@@ -2425,18 +2425,18 @@ static int vcf_parse_format(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v, char *p
             } else if ((z->y>>4&0xf) == BCF_HT_INT) {
                 int32_t *x = (int32_t*)(z->buf + z->size * (size_t)m);
                 for (l = 0;; ++t) {
-                    if (*t == '.') x[l++] = bcf_int32_missing, ++t; // ++t to skip "."
-                    else
-                    {
+                    if (*t == '.') {
+                        x[l++] = bcf_int32_missing, ++t; // ++t to skip "."
+                    } else {
                         overflow = 0;
                         char *te;
                         long int tmp_val = hts_str2int(t, &te, sizeof(tmp_val)*CHAR_BIT, &overflow);
                         if ( te==t || overflow || tmp_val<BCF_MIN_BT_INT32 || tmp_val>BCF_MAX_BT_INT32 )
                         {
-                            if ( !extreme_int_warned )
+                            if ( !extreme_val_warned )
                             {
-                                hts_log_warning("Extreme FORMAT/%s value encountered and set to missing at %s:%"PRIhts_pos,h->id[BCF_DT_ID][fmt[j-1].key].key,bcf_seqname_safe(h,v), v->pos+1);
-                                extreme_int_warned = 1;
+                                hts_log_warning("Extreme FORMAT/%s value encountered and set to missing at %s:%"PRIhts_pos, h->id[BCF_DT_ID][fmt[j-1].key].key, bcf_seqname_safe(h,v), v->pos+1);
+                                extreme_val_warned = 1;
                             }
                             tmp_val = bcf_int32_missing;
                         }
@@ -2450,8 +2450,20 @@ static int vcf_parse_format(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v, char *p
             } else if ((z->y>>4&0xf) == BCF_HT_REAL) {
                 float *x = (float*)(z->buf + z->size * (size_t)m);
                 for (l = 0;; ++t) {
-                    if (*t == '.' && !isdigit_c(t[1])) bcf_float_set_missing(x[l++]), ++t; // ++t to skip "."
-                    else x[l++] = strtod(t, &t);
+                    if (*t == '.' && !isdigit_c(t[1])) {
+                        bcf_float_set_missing(x[l++]), ++t; // ++t to skip "."
+                    } else {
+                        overflow = 0;
+                        char *te;
+                        float tmp_val = hts_str2dbl(t, &te, &overflow);
+                        if ( (te==t || overflow) && !extreme_val_warned )
+                        {
+                            hts_log_warning("Extreme FORMAT/%s value encountered at %s:%"PRIhts_pos, h->id[BCF_DT_ID][fmt[j-1].key].key, bcf_seqname(h,v), v->pos+1);
+                            extreme_val_warned = 1;
+                        }
+                        x[l++] = tmp_val;
+                        t = te;
+                    }
                     if (*t != ',') break;
                 }
                 if ( !l ) bcf_float_set_missing(x[l++]);    // An empty field, insert missing value
@@ -2768,8 +2780,9 @@ static int vcf_parse_info(kstring_t *str, const bcf_hdr_t *h, bcf1_t *v, char *p
                 float *val_f = (float *)a_val;
                 for (i = 0, t = val; i < n_val; ++i, ++t)
                 {
-                    val_f[i] = strtod(t, &te);
-                    if ( te==t ) // conversion failed
+                    overflow = 0;
+                    val_f[i] = hts_str2dbl(t, &te, &overflow);
+                    if ( te==t || overflow ) // conversion failed
                         bcf_float_set_missing(val_f[i]);
                     for (t = te; *t && *t != ','; t++);
                 }