Update get_blocks function

src/branch_sequences.c

@@ -744,6 +744,8 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i
 
     // create the pileup of the snps and their sphere of influence
     int snp_counter = 0;
+    int min_postion = genome_size;
+    int max_position = 0;
     for(snp_counter = 0; snp_counter < number_of_branch_snps; snp_counter++)
     {
         int j = 0;
@@ -760,6 +762,11 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i
             snp_sliding_window_counter = 0;
         }
 
+        if(snp_sliding_window_counter < min_postion)
+        {
+            min_postion = snp_sliding_window_counter;
+        }
+
         // Upper bound of the window around a snp
         int max_snp_sliding_window_counter = snp_site_coords[snp_counter]+(window_size/2);
         max_snp_sliding_window_counter = extend_upper_part_of_window(snp_site_coords[snp_counter] + 1,
@@ -771,6 +778,11 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i
         {
             max_snp_sliding_window_counter = genome_size;
         }
+
+        if(max_snp_sliding_window_counter > max_position)
+        {
+          max_position= max_snp_sliding_window_counter;
+        }
 
         for(j = snp_sliding_window_counter; j < max_snp_sliding_window_counter; j++)
         {
@@ -784,7 +796,7 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i
     int block_lower_bound = 0;
     // Scan across the pileup and record where blocks are above the cutoff
     int i;
-    for(i = 0; i <= genome_size; i++)
+    for(i = min_postion; i <= max_position; i++)
     {
         // Just entered the start of a block
         if(window_count[i] > cutoff && in_block == 0)
@@ -794,20 +806,23 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i
         }
 
         // Reached end of genome
-        if(i == genome_size && in_block == 1)
-        {
-            block_coordinates[0][number_of_blocks] = block_lower_bound;
-            block_coordinates[1][number_of_blocks] = i;
-            number_of_blocks++;
-            in_block = 0;
-        }
-        // Just left a block
-        else if(window_count[i] <= cutoff && in_block == 1)
+        if (in_block == 1)
         {
-            block_coordinates[0][number_of_blocks] = block_lower_bound;
-            block_coordinates[1][number_of_blocks] = i-1;
-            number_of_blocks++;
-            in_block = 0;
+          if(i == genome_size)
+          {
+              block_coordinates[0][number_of_blocks] = block_lower_bound;
+              block_coordinates[1][number_of_blocks] = i;
+              number_of_blocks++;
+              in_block = 0;
+          }
+          // Just left a block
+          else if(window_count[i] <= cutoff)
+          {
+              block_coordinates[0][number_of_blocks] = block_lower_bound;
+              block_coordinates[1][number_of_blocks] = i-1;
+              number_of_blocks++;
+              in_block = 0;
+          }
         }
 
     }