update to work with kmer rarity

1 files changed, 120 insertions, 11 deletions

diff --git a/src/c/quikr.c b/src/c/quikr.c
index 9fc0759..1b8daa3 100644
--- a/src/c/quikr.c
+++ b/src/c/quikr.c
@@ -1,4 +1,3 @@
-#include <ctype.h>
 #include <errno.h>
 #include <getopt.h>
 #include <math.h>
@@ -24,6 +23,14 @@ int main(int argc, char **argv) {
   char *output_filename = NULL;
 
   unsigned long long x = 0;
+  unsigned long long y = 0;
+  unsigned long long z = 0;
+
+  unsigned long long rare_value = 0;
+  unsigned long long rare_width = 0; 
+
+  double percentage = 1.0;
+  double rare_percent = 1.0;
 
   unsigned long long width = 0;
 
@@ -39,6 +46,7 @@ int main(int argc, char **argv) {
       {"lambda",  required_argument, 0, 'l'},
       {"output", required_argument, 0, 'o'},
       {"sensing-matrix", required_argument, 0, 's'},
+      {"rare-percent", required_argument, 0, 'r'},
       {"verbose", no_argument, 0, 'v'},
       {"version", no_argument, 0, 'V'},
       {"help", no_argument, 0, 'h'},
@@ -48,7 +56,7 @@ int main(int argc, char **argv) {
 
     int option_index = 0;
 
-    c = getopt_long (argc, argv, "k:l:s:i:o:hdvV", long_options, &option_index);
+    c = getopt_long (argc, argv, "k:l:s:r:i:o:r:hdvV", long_options, &option_index);
 
     if (c == -1)
       break;
@@ -60,6 +68,9 @@ int main(int argc, char **argv) {
       case 'l':
         lambda = atoi(optarg);
         break;
+      case 'r':
+        rare_percent = atof(optarg);
+        break;
       case 's':
         sensing_matrix_filename = optarg;
         break;
@@ -75,17 +86,14 @@ int main(int argc, char **argv) {
       case 'V':
         printf("%s\n", VERSION);
         exit(EXIT_SUCCESS);
-        break;
       case 'h':
         printf("%s\n", USAGE);
         exit(EXIT_SUCCESS);
-        break;
       default:
         break;
     }
   }
 
-
   if(sensing_matrix_filename == NULL) {
     fprintf(stderr, "Error: sensing matrix filename (-s) must be specified\n\n");
     fprintf(stderr, "%s\n", USAGE);
@@ -102,7 +110,8 @@ int main(int argc, char **argv) {
     exit(EXIT_FAILURE);
   }
 
-  if(verbose) { 
+  if(verbose) {
+    printf("rare width:%ld\n", rare_width);
     printf("kmer: %u\n", kmer);
     printf("lambda: %llu\n", lambda);
     printf("fasta: %s\n", input_fasta_filename);
@@ -124,18 +133,109 @@ int main(int argc, char **argv) {
     exit(EXIT_FAILURE);
 	}
 
+
   // 4 "ACGT" ^ Kmer gives us the size of output rows
-  width = pow(4, kmer);
-  width = width + 1;
+  width = pow_four(kmer);
 
 	// load counts matrix
-	double *count_matrix = setup_count_matrix(input_fasta_filename, kmer, lambda, width); 
+	double *count_matrix = malloc(width * sizeof(double));
+	check_malloc(count_matrix, NULL);
+
+	{
+		unsigned long long *integer_counts = get_kmer_counts_from_file(input_fasta_filename, kmer);
+
+		for(x = 0; x < width; x++)
+			count_matrix[x] = (double)integer_counts[x];
+
+		free(integer_counts);
+	}
 
 	// load sensing matrix
-	struct matrix *sensing_matrix = setup_sensing_matrix(sensing_matrix_filename, kmer, lambda, width); 
+	struct matrix *sensing_matrix = load_sensing_matrix(sensing_matrix_filename, kmer);
+
+
+	// get our "rare" counts
+	while(1) {
+		rare_width = 0;
+		for(x = 0; x < width; x++) {
+			if(count_matrix[x] < rare_value) {
+				rare_width++;
+			}
+		}
+		percentage = (float)rare_width / (float)width;
+
+		if(percentage >= rare_percent)
+			break;
+
+		rare_value++;
+	}
+
+	if(verbose)
+		printf("there are %llu values less than %llu\n", rare_width, rare_value);
+
+	// add a extra space for our zero's array
+	rare_width++;
+
+	// store our count matrix
+	double *count_matrix_rare = calloc(rare_width, sizeof(double));
+	check_malloc(count_matrix_rare, NULL);
+
+	double *sensing_matrix_rare = malloc((rare_width) * sensing_matrix->sequences * sizeof(double));
+	check_malloc(sensing_matrix_rare, NULL);
+
+	// copy only kmers from our original counts that match our rareness percentage
+	// in both our count matrix and our sensing matrix
+	for(x = 0, y = 1;  x < width; x++) {
+		if(count_matrix[x] < rare_value) {
+			count_matrix_rare[y] = count_matrix[x];
+
+			for(z = 0; z < sensing_matrix->sequences; z++) 
+				sensing_matrix_rare[z*rare_width + y] = sensing_matrix->matrix[z*width + x];
+
+			y++;
+		}
+	}
+
+	// normalize our kmer counts and our sensing_matrix
+	normalize_matrix(count_matrix_rare, 1, rare_width);
+	normalize_matrix(sensing_matrix_rare, sensing_matrix->sequences, rare_width);
+
+	// multiply our kmer counts and sensing matrix by lambda
+	for(x = 1; x < rare_width; x++) 
+		count_matrix_rare[x] *= lambda;
+
+	for(x = 0; x < sensing_matrix->sequences; x++) {
+		for(y = 1; y < rare_width; y++) {
+			sensing_matrix_rare[rare_width*x + y] *= lambda;
+		}
+	}
+
+	// count_matrix's first element should be zero
+	count_matrix_rare[0] = 0;
+	// stack one's on our first row of our sensing matrix
+	for(x = 0; x < sensing_matrix->sequences; x++) {
+		sensing_matrix_rare[x*rare_width] = 1.0;
+	}
+
+	// DEBUG OUR MATRIX
+	if(verbose) {
+		printf("COUNT_MATRIX\n");
+		for(x = 0; x < rare_width; x++) 
+			printf("%lf\t", count_matrix_rare[x]);
+			printf("\n");
+
+		printf("SENSING_MATRIX\n");
+		for(y = 0; y < sensing_matrix->sequences - 1; y++) {
+			for(x = 0; x < rare_width; x++) {
+				printf("%lf\t", sensing_matrix_rare[x*rare_width + y]);
+			}
+			printf("%lf\n", sensing_matrix_rare[x*rare_width + (rare_width - 1)]);
+		}
+
+	}
 
 	// run NNLS
-  double *solution = nnls(sensing_matrix->matrix, count_matrix, sensing_matrix->sequences, width);
+  double *solution = nnls(sensing_matrix_rare, count_matrix_rare, sensing_matrix->sequences, rare_width);
 
   // normalize our solution vector
   normalize_matrix(solution, 1, sensing_matrix->sequences);
@@ -149,7 +249,16 @@ int main(int argc, char **argv) {
   for(x = 0; x < sensing_matrix->sequences; x++) {
       fprintf(output_fh, "%.10lf\n", solution[x]);
   }
+
   fclose(output_fh);
 
+	free(solution);
+
+	free(count_matrix);
+	free(sensing_matrix);
+
+	free(count_matrix_rare);
+	free(sensing_matrix_rare);
+
   return EXIT_SUCCESS;
 }