-
Notifications
You must be signed in to change notification settings - Fork 12
/
compareDistributions.c
107 lines (80 loc) · 3.51 KB
/
compareDistributions.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
// Compare nanopore HDP distributions
#include <stdio.h>
#include <inttypes.h>
#include <stdbool.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <string.h>
#include "sonLib.h"
#include "nanopore_hdp.h"
#include "hdp_math_utils.h"
#define FILENAME_BUFFER_LEN 100
char* kmer_from_index(int64_t index, const char* alphabet, int64_t alphabet_size, int64_t kmer_length) {
char* kmer = (char*) malloc(sizeof(char) * (kmer_length + 1));
int64_t index_remainder = index;
kmer[kmer_length] = '\0';
for (int64_t i = kmer_length - 1; i >= 0; i--) {
kmer[i] = alphabet[index_remainder % alphabet_size];
index_remainder /= alphabet_size;
}
return kmer;
}
void write_kmer_distr(NanoporeHDP* nhdp, char* kmer, double* eval_grid, int64_t grid_length, char *workingDirectory) {
//char filename[FILENAME_BUFFER_LEN];
char *filename = stString_print("%s/%s_distr.txt", workingDirectory, kmer);
//sprintf(filename, "%s/%s_distr.txt", workingDirectory, kmer);
FILE* out = fopen(filename, "w");
//fprintf(out, "%s\t", kmer);
double density;
for (int64_t i = 0; i < grid_length - 1; i++) {
density = get_nanopore_kmer_density(nhdp, kmer, eval_grid + i);
fprintf(out, "%.17lg\n", density);
}
density = get_nanopore_kmer_density(nhdp, kmer, eval_grid + (grid_length - 1));
fprintf(out, "%.17lg\n", density);
fclose(out);
}
void write_all_kmer_distrs(NanoporeHDP* nhdp, double* eval_grid, int64_t grid_length, char *workingDirectory) {
char *x_filename = stString_print("%s/x_vals.txt", workingDirectory);
//char *distributions = stString_print("%s/distributions.tsv", workingDirectory);
//sprintf(x_filename, "%s/x_vals.txt", workingDirectory);
//sprintf(distributions, "%s/distributions.tsv", workingDirectory);
FILE* x_vals = fopen(x_filename, "w");
//FILE *distributionFileHandle = fopen(distributions, "a");
for (int64_t i = 0; i < grid_length - 1; i++) {
fprintf(x_vals, "%.17lg\n", eval_grid[i]);
}
fprintf(x_vals, "%.17lg", eval_grid[grid_length - 1]);
fclose(x_vals);
int64_t alphabet_size = get_nanopore_hdp_alphabet_size(nhdp);
int64_t kmer_length = get_nanopore_hdp_kmer_length(nhdp);
char* alphabet = get_nanopore_hdp_alphabet(nhdp);
int64_t num_kmers = power(alphabet_size, kmer_length);
for (int64_t kmer_index = 0; kmer_index < num_kmers; kmer_index++) {
char* kmer = kmer_from_index(kmer_index, alphabet, alphabet_size, kmer_length);
write_kmer_distr(nhdp, kmer, eval_grid, grid_length, workingDirectory);
free(kmer);
}
free(alphabet);
//fclose(workingDirectory);
}
int main(int argc, char *argv[]) {
if (argc != 3) {
fprintf(stderr, "USAGE_NEW: compareDistributions [NanoporeHDP_file] [distribution_directory]\n");
exit(EXIT_FAILURE);
}
char *modelFile = argv[1];
fprintf(stderr, "[compareDistributions] NOTICE: Loading NanoporeHDP from %s\n", modelFile);
char *workingDirectory = argv[2];
//FILE *fH = fopen(destinationFile, "w");
fprintf(stderr, "[compareDistributions] NOTICE: Putting distributions in %s\n", workingDirectory);
NanoporeHDP *nHdp= deserialize_nhdp(modelFile);
double gridStart = 30.0;
double gridStop = 90.0;
int64_t gridLength = 600;
double *evaluationGrid = linspace(gridStart, gridStop, gridLength);
write_all_kmer_distrs(nHdp, evaluationGrid, gridLength, workingDirectory);
free(evaluationGrid);
destroy_nanopore_hdp(nHdp);
return 0;
}