Improve comments and self-describing code

Author: jkbonfield

vcf.c

@@ -2663,6 +2663,25 @@ typedef struct {
     uint8_t *buf;       // Pointer into h->mem
 } fmt_aux_t;
 
+// fmt_aux_t field notes:
+// max_* are biggest sizes of the various FORMAT fields across all samples.
+// We use these after pivoting the data to ensure easy random access
+// of a specific sample.
+//
+// max_m is only used for type BCF_HT_REAL or BCF_HT_INT
+// max_g is only used for is_gt == 1 (will be BCF_HT_STR)
+// max_l is only used for is_gt == 0 (will be BCF_HT_STR)
+//
+// These are computed in vcf_parse_format_max3 and used in
+// vcf_parse_format_alloc4 to get the size.
+//
+// size is computed from max_g, max_l, max_m and is_gt.  Once computed
+// the max values are never accessed again.
+//
+// In theory all 4 vars could be coalesced into a single variable, but this
+// significantly harms speed (even if done via a union).  It's about 25-30%
+// slower.
+
 static inline int align_mem(kstring_t *s)
 {
     int e = 0;
@@ -2919,15 +2938,19 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
         while ( t < end )
         {
             fmt_aux_t *z = &fmt[j++];
+            const int htype = z->y>>4&0xf;
             if (!z->buf) {
                 hts_log_error("Memory allocation failure for FORMAT field type %d at %s:%"PRIhts_pos,
                               z->y>>4&0xf, bcf_seqname_safe(h,v), v->pos+1);
                 v->errcode |= BCF_ERR_LIMITS;
                 return -1;
             }
-            if ((z->y>>4&0xf) == BCF_HT_STR) {
+
+            if (htype == BCF_HT_STR) {
                 int l;
-                if (z->is_gt) { // genotypes
+                if (z->is_gt) {
+                    // Genotypes.
+                    // <val>([|/]<val>)+... where <val> is [0-9]+ or ".".
                     int32_t is_phased = 0;
                     uint32_t *x = (uint32_t*)(z->buf + z->size * (size_t)m);
                     uint32_t unreadable = 0;
@@ -2966,14 +2989,20 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
                         return -1;
                     }
                     if ( !l ) x[l++] = 0;   // An empty field, insert missing value
-                    for (; l < z->size>>2; ++l) x[l] = bcf_int32_vector_end;
+                    for (; l < z->size>>2; ++l)
+                        x[l] = bcf_int32_vector_end;
+
                 } else {
+                    // Otherwise arbitrary strings
                     char *x = (char*)z->buf + z->size * (size_t)m;
-                    for (l = 0; *t != ':' && *t; ++t) x[l++] = *t;
+                    for (l = 0; *t != ':' && *t; ++t)
+                        x[l++] = *t;
                     if (z->size > l)
                         memset(&x[l], 0, (z->size-l) * sizeof(*x));
                 }
-            } else if ((z->y>>4&0xf) == BCF_HT_INT) {
+
+            } else if (htype == BCF_HT_INT) {
+                // One or more integers in an array
                 int32_t *x = (int32_t*)(z->buf + z->size * (size_t)m);
                 int l;
                 for (l = 0;; ++t) {
@@ -2987,7 +3016,8 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
                         {
                             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);
+                                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;
@@ -2997,9 +3027,13 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
                     }
                     if (*t != ',') break;
                 }
-                if ( !l ) x[l++] = bcf_int32_missing;
-                for (; l < z->size>>2; ++l) x[l] = bcf_int32_vector_end;
-            } else if ((z->y>>4&0xf) == BCF_HT_REAL) {
+                if ( !l )
+                    x[l++] = bcf_int32_missing;
+                for (; l < z->size>>2; ++l)
+                    x[l] = bcf_int32_vector_end;
+
+            } else if (htype == BCF_HT_REAL) {
+                // One of more floating point values in an array
                 float *x = (float*)(z->buf + z->size * (size_t)m);
                 int l;
                 for (l = 0;; ++t) {
@@ -3019,10 +3053,13 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
                     }
                     if (*t != ',') break;
                 }
-                if ( !l ) bcf_float_set_missing(x[l++]);    // An empty field, insert missing value
-                for (; l < z->size>>2; ++l) bcf_float_set_vector_end(x[l]);
+                if ( !l )
+                    // An empty field, insert missing value
+                    bcf_float_set_missing(x[l++]);
+                for (; l < z->size>>2; ++l)
+                    bcf_float_set_vector_end(x[l]);
             } else {
-                hts_log_error("Unknown FORMAT field type %d at %s:%"PRIhts_pos, z->y>>4&0xf, bcf_seqname_safe(h,v), v->pos+1);
+                hts_log_error("Unknown FORMAT field type %d at %s:%"PRIhts_pos, htype, bcf_seqname_safe(h,v), v->pos+1);
                 v->errcode |= BCF_ERR_TAG_INVALID;
                 return -1;
             }
@@ -3043,10 +3080,12 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
             }
         }
 
-        for (; j < v->n_fmt; ++j) { // fill end-of-vector values
+        // fill end-of-vector values
+        for (; j < v->n_fmt; ++j) {
             fmt_aux_t *z = &fmt[j];
+            const int htype = z->y>>4&0xf;
             int l;
-            if ((z->y>>4&0xf) == BCF_HT_STR) {
+            if (htype == BCF_HT_STR) {
                 if (z->is_gt) {
                     int32_t *x = (int32_t*)(z->buf + z->size * (size_t)m);
                     if (z->size) x[0] = bcf_int32_missing;
@@ -3058,11 +3097,11 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
                         memset(&x[1], 0, (z->size-1) * sizeof(*x));
                     }
                 }
-            } else if ((z->y>>4&0xf) == BCF_HT_INT) {
+            } else if (htype == BCF_HT_INT) {
                 int32_t *x = (int32_t*)(z->buf + z->size * (size_t)m);
                 x[0] = bcf_int32_missing;
                 for (l = 1; l < z->size>>2; ++l) x[l] = bcf_int32_vector_end;
-            } else if ((z->y>>4&0xf) == BCF_HT_REAL) {
+            } else if (htype == BCF_HT_REAL) {
                 float *x = (float*)(z->buf + z->size * (size_t)m);
                 bcf_float_set_missing(x[0]);
                 for (l = 1; l < z->size>>2; ++l) bcf_float_set_vector_end(x[l]);
@@ -3145,6 +3184,23 @@ static int vcf_parse_format(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
     if (vcf_parse_format_dict2(s, h, v, p, q, fmt) < 0)
         return -1;
 
+    // FORMAT data is per-sample A:B:C A:B:C A:B:C ... but in memory it is
+    // stored as per-type arrays AAA... BBB... CCC...  This is basically
+    // a data rotation or pivot.
+
+    // The size of elements in the array grow to their maximum needed,
+    // permitting fast random access.  This means however we have to first
+    // scan the whole FORMAT line to find the maximum of each type, and
+    // then scan it again to find the store the data.
+    // We break this down into compute-max, allocate, fill-out-buffers
+
+    // TODO: ?
+    // The alternative would be to pivot on the first pass, with fixed
+    // size entries for numerics and concatenated strings otherwise, also
+    // tracking maximum sizes.  Then on a second pass we reallocate and
+    // copy the data again to a uniformly sized array.  Two passes through
+    // memory, but without doubling string parsing.
+
     // compute max
     if (vcf_parse_format_max3(s, h, v, p, q, fmt) < 0)
         return -1;