tests are passing with speed-up improvements, a few improvements are left out to make test pass for now

Author: eneskuluk

htslib/vcf.h

@@ -1605,8 +1605,9 @@ static inline int bcf_enc_size(kstring_t *s, size_t size, int type)
             s->l += 6;
         }
         else{
-            *p++ = 1<<4|BCF_BT_INT64;
-            s->l += 10;
+            *p++ = 1<<4|BCF_BT_INT64;//
+            i64_to_le(size,p);
+            s->l += 10; // not so sure about +10 here, whether it is accurate or changes anything.
             return -1;
         }
     }

synced_bcf_reader.c

@@ -1048,7 +1048,6 @@ static bcf_sr_regions_t *_regions_init_string(const char *str)
     kstring_t tmp = {0,0,0};
     const char *sp = str, *ep = str;
     hts_pos_t from, to;
-    unsigned char inside_quotes = 0;
     while ( 1 )
     {
         tmp.l = 0;
@@ -1065,28 +1064,8 @@ static bcf_sr_regions_t *_regions_init_string(const char *str)
         }
         else
         {
-        //A quote is seen, flip flag inside_quotes
-        if(*ep == '"')
-        {
-            inside_quotes = 1 ^ inside_quotes;
-            sp = ++ep;
-        }
-        while ( *ep && ((inside_quotes && *ep!='"') || (!inside_quotes && *ep!=',' && *ep!=':')) ) ep++;
-        tmp.l = 0;
-        kputsn(sp,ep-sp,&tmp);
-        if(inside_quotes)
-        {
-            if(*ep == '"')
-            {
-                inside_quotes = 0;
-                ++ep;
-            }
-            else
-            {
-                fprintf(stderr,"[%s:%d %s] Could not parse the region(s): %s - terminating \" missing\n", __FILE__,__LINE__,__FUNCTION__,str);
-                free(reg); free(tmp.s); return NULL;
-            }
-        }
+            while ( *ep && *ep!=',' && *ep!=':' ) ep++;
+            kputsn(sp,ep-sp,&tmp);
         }
         if ( *ep==':' )
         {

vcf.c

@@ -2621,29 +2621,68 @@ static int bcf_enc_long1(kstring_t *s, int64_t x) {
 }
 #endif
 
-int bcf_enc_vlong(kstring_t *s, const int n, const int64_t *a, int wsize)
+int bcf_enc_vlong(kstring_t *s, const int n, const int64_t *a, int wsize)//supposed to be optimized, not working yet.
 {
     int64_t max = INT64_MIN, min = INT64_MAX;
     int i;
-    if (n <= 0) return bcf_enc_size(s, 0, BCF_BT_NULL);
-    else if (n == 1) return bcf_enc_long1(s, a[0]);
-    else {
+    if (n <= 0) {
+        return bcf_enc_size(s, 0, BCF_BT_NULL);
+    } else if (n == 1) {
+        return bcf_enc_long1(s, a[0]);
+    } else {
         if (wsize <= 0) wsize = n;
-        for (i = 0; i < n; ++i) {
-            if (a[i] == bcf_int64_missing || a[i] == bcf_int64_vector_end ) continue;
+
+        // Equivalent to:
+        // for (i = 0; i < n; ++i) {
+        //     if (a[i] == bcf_int32_missing || a[i] == bcf_int32_vector_end )
+        //         continue;
+        //     if (max < a[i]) max = a[i];
+        //     if (min > a[i]) min = a[i];
+        // }
+        int64_t max4[4] = {INT64_MIN, INT64_MIN, INT64_MIN, INT64_MIN};
+        int64_t min4[4] = {INT64_MAX, INT64_MAX, INT64_MAX, INT64_MAX};
+        for (i = 0; i < (n&~3); i+=4) {
+            // bcf_int32_missing    == INT32_MIN and
+            // bcf_int32_vector_end == INT32_MIN+1.
+            // We skip these, but can mostly avoid explicit checking
+            if (max4[0] < a[i+0]) max4[0] = a[i+0];
+            if (max4[1] < a[i+1]) max4[1] = a[i+1];
+            if (max4[2] < a[i+2]) max4[2] = a[i+2];
+            if (max4[3] < a[i+3]) max4[3] = a[i+3];
+            if (min4[0] > a[i+0] && a[i+0] > INT64_MIN+1) min4[0] = a[i+0];
+            if (min4[1] > a[i+1] && a[i+1] > INT64_MIN+1) min4[1] = a[i+1];
+            if (min4[2] > a[i+2] && a[i+2] > INT64_MIN+1) min4[2] = a[i+2];
+            if (min4[3] > a[i+3] && a[i+3] > INT64_MIN+1) min4[3] = a[i+3];
+        }
+        min = min4[0];
+        if (min > min4[1]) min = min4[1];
+        if (min > min4[2]) min = min4[2];
+        if (min > min4[3]) min = min4[3];
+        max = max4[0];
+        if (max < max4[1]) max = max4[1];
+        if (max < max4[2]) max = max4[2];
+        if (max < max4[3]) max = max4[3];
+        for (; i < n; ++i) {
             if (max < a[i]) max = a[i];
-            if (min > a[i]) min = a[i];
+            if (min > a[i] && a[i] > INT64_MIN+1) min = a[i];
         }
+
         if (max <= BCF_MAX_BT_INT8 && min >= BCF_MIN_BT_INT8) {
-            bcf_enc_size(s, wsize, BCF_BT_INT8);
-            for (i = 0; i < n; ++i)
-                if ( a[i]==bcf_int64_vector_end ) kputc(bcf_int8_vector_end, s);
-                else if ( a[i]==bcf_int64_missing ) kputc(bcf_int8_missing, s);
-                else kputc(a[i], s);
+            if (bcf_enc_size(s, wsize, BCF_BT_INT8) < 0 ||
+                ks_resize(s, s->l + n) < 0)
+                return -1;
+            uint8_t *p = (uint8_t *) s->s + s->l;
+            for (i = 0; i < n; ++i, p++) {
+                if ( a[i]==bcf_int64_vector_end )   *p = bcf_int8_vector_end;
+                else if ( a[i]==bcf_int64_missing ) *p = bcf_int8_missing;
+                else *p = a[i];
+            }
+            s->l += n;
         } else if (max <= BCF_MAX_BT_INT16 && min >= BCF_MIN_BT_INT16) {
             uint8_t *p;
-            bcf_enc_size(s, wsize, BCF_BT_INT16);
-            ks_resize(s, s->l + n * sizeof(int16_t));
+            if (bcf_enc_size(s, wsize, BCF_BT_INT16) < 0 ||
+                ks_resize(s, s->l + n * sizeof(int16_t)) < 0)
+                return -1;
             p = (uint8_t *) s->s + s->l;
             for (i = 0; i < n; ++i)
             {
@@ -2657,8 +2696,9 @@ int bcf_enc_vlong(kstring_t *s, const int n, const int64_t *a, int wsize)
             s->l += n * sizeof(int16_t);
         } else if(max <= BCF_MAX_BT_INT32 && min >= BCF_MIN_BT_INT32){
             uint8_t *p;
-            bcf_enc_size(s, wsize, BCF_BT_INT32);
-            ks_resize(s, s->l + n * sizeof(int32_t));
+            if (bcf_enc_size(s, wsize, BCF_BT_INT32) < 0 ||
+                ks_resize(s, s->l + n * sizeof(int32_t)) < 0)
+                return -1;
             p = (uint8_t *) s->s + s->l;
             for (i = 0; i < n; ++i) {
 		int32_t x;
@@ -2670,11 +2710,11 @@ int bcf_enc_vlong(kstring_t *s, const int n, const int64_t *a, int wsize)
             }
             s->l += n * sizeof(int32_t);
         }
-#ifdef VCF_ALLOW_INT64
+        #ifdef VCF_ALLOW_INT64
         else {
             uint8_t *p;
-            bcf_enc_size(s, wsize, BCF_BT_INT64);
-            ks_resize(s, s->l + n * sizeof(int64_t));
+            if(bcf_enc_size(s, wsize, BCF_BT_INT64) < 0 || ks_resize(s, s->l + n * sizeof(int64_t)) < 0)
+            return -1;
             p = (uint8_t *) s->s + s->l;
             for (i = 0; i < n; ++i) {
 		int64_t x = a[i];
@@ -2688,7 +2728,7 @@ int bcf_enc_vlong(kstring_t *s, const int n, const int64_t *a, int wsize)
 #endif
     }
 
-    return 0; // FIXME: check for errs in this function
+    return 0;
 }
 
 static inline int serialize_float_array(kstring_t *s, size_t n, const float *a) {
@@ -3524,7 +3564,7 @@ static int vcf_parse_info(kstring_t *str, const bcf_hdr_t *h, bcf1_t *v, char *p
                 v->unpacked |= BCF_IS_64BIT;
                 bcf_enc_vlong(str, n_val, a_val, -1);
                 val1 = a_val[0];
-                if (n_val==1 && strcmp(key, "END") == 0)
+                if (n_val==1 && strcmp(key, "END") == 0)//memset instead of strcmp
                 {
                     if ( val1 <= v->pos )
                     {
@@ -3586,6 +3626,7 @@ int vcf_parse(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v)
     // Ensure string we parse has space to permit some over-flow when during
     // parsing.  Eg to do memcmp(key, "END", 4) in vcf_parse_info over
     // the more straight forward looking strcmp, giving a speed advantage.
+    /*
     if (ks_resize(s, s->l+4) < 0)
         return -1;
 
@@ -3598,7 +3639,7 @@ int vcf_parse(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v)
     s->s[s->l+1] = 0;
     s->s[s->l+2] = 0;
     s->s[s->l+3] = 0;
-
+    */ // commented out the part that was required for optimization in vcf_parse_info function, will take a look later.
     bcf_clear1(v);
     str = &v->shared;
     memset(&aux, 0, sizeof(ks_tokaux_t));
@@ -4657,7 +4698,7 @@ int bcf_hdr_parse_required_sample_line(bcf_hdr_t *hdr, char *htxt, size_t* hdr_l
         // operations do not really care about a few malformed lines).
         // In the future we may want to add a strict mode that errors in
         // this case.
-        if ( strncmp("#CHROM\tPOS",p,10) != 0 ) {
+        if ( strncmp("#CHROM\t",p,7) && strncmp("#CHROM ",p,7) ) {
             char *eol = strchr(p, '\n');
             if (*p != '\0') {
                 char buffer[320];