Revert breaking changes

src/branch_sequences.c

@@ -77,27 +77,33 @@ int copy_and_concat_2d_integer_arrays(int ** array_1, int array_1_size, int ** a
 // Loop over all recombination blocks and return a list of all snp indices that fall within those blocks.
 int get_list_of_snp_indices_which_fall_in_downstream_recombinations(int ** current_block_coordinates,int num_blocks, int * snp_locations,int number_of_snps, int * snps_in_recombinations)
 {
-  int num_snps_in_recombinations =0;
-  int i = 0;
+	int num_snps_in_recombinations =0;
+	int i = 0;
 
   // loop over each block
-  for(i = 0; i<num_blocks; i++ )
-  {
-    int current_index = 0;
+	for(i = 0; i<num_blocks; i++ )
+	{
+		int current_index = 0;
     // convert the starting coordinates of block to the nearest SNP index
-    current_index = find_starting_index(current_block_coordinates[0][i],snp_locations,0, number_of_snps);
-    
-    // starting at the begining index of block, count all the snps until the end of the block.
-    int j;
-    for(j = current_index; (j < number_of_snps && snp_locations[j] <= current_block_coordinates[1][i]); j++)
+		current_index = find_starting_index(current_block_coordinates[0][i],snp_locations,0, number_of_snps);
+		
+    //make sure that the index begins at start of block
+		int beginning_j = current_index;
+    for(beginning_j = current_index; snp_locations[beginning_j] < current_block_coordinates[0][i];beginning_j++)
     {
-        snps_in_recombinations[num_snps_in_recombinations] = j;
-        num_snps_in_recombinations++;
     }
-  }
+
+    int j;
+    // starting at the begining index of block, count all the snps until the end of the bock.
+		for(j = beginning_j; (j < number_of_snps && snp_locations[j] <= current_block_coordinates[1][i]); j++)
+		{
+				snps_in_recombinations[num_snps_in_recombinations] = j;
+				num_snps_in_recombinations++;
+		}
+	}
 
   // may contain duplications
-  return num_snps_in_recombinations;
+	return num_snps_in_recombinations;
 }
 
 
@@ -738,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;
@@ -754,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,
@@ -765,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++)
         {
@@ -778,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)
@@ -788,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;
+          }
         }
 
     }