+ added small app to ease profiling
[qpalma.git] / dyn_prog / test_fill_matrix.cpp
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <fstream>
4 #include <sstream>
5 using namespace std;
6
7 #include "qpalma_dp.h"
8
9 void initialize_quality_plifs(penalty_struct* qplifs, double* params) {
10
11 printf("begin of initialize_quality_plifs \n");
12 double limits[10] = {-5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0};
13
14 int idx, pos;
15 int ctr = 36;
16
17 for(idx=0;idx<30;idx++) {
18 penalty_struct new_plif;
19 init_penalty_struct(new_plif);
20
21 new_plif.min_len = -5;
22 new_plif.max_len = 40;
23 new_plif.len = 10;
24
25 new_plif.limits = (double*) malloc(sizeof(double)*10);
26 if (new_plif.limits == NULL)
27 printf("malloc (limits)\n");
28
29 new_plif.penalties = (double*) malloc(sizeof(double)*10);
30 if (new_plif.penalties == NULL)
31 printf("malloc (penalties)\n");
32
33 for(pos=0;pos<10;pos++) {
34 new_plif.limits[pos] = limits[pos];
35 new_plif.penalties[pos] = params[ctr];
36 ctr++;
37 }
38
39 qplifs[idx] = new_plif;
40 }
41
42 printf("end of initialize_quality_plifs \n");
43 return;
44 }
45
46
47 int main(int argc, char** argv){
48
49 printf("Starting test_fm\n");
50 int idx;
51 char* alphabet = (char*) malloc(sizeof(char)*4);
52
53 alphabet[0] = 'a';
54 alphabet[1] = 'c';
55 alphabet[2] = 'g';
56 alphabet[3] = 't';
57
58 int est_len = 36;
59 int dna_len = 10000;
60 int nr_paths = 2;
61
62 int nr_donor_sites = 0 ;
63 int d_len = dna_len;
64
65 printf("begin of donor/acceptor stuff\n");
66 double* donor = (double*) malloc(sizeof(double)*dna_len);
67 double* acceptor = (double*) malloc(sizeof(double)*dna_len);
68
69 for(idx=0;idx<dna_len;idx++) {
70 donor[idx] = drand48();
71 acceptor[idx] = drand48();
72 }
73
74 for (int ii=0; ii<d_len; ii++) {
75 if(isfinite(donor[ii])) {
76 nr_donor_sites++ ;
77 }
78 }
79
80 int* donor_sites = new int[nr_donor_sites];
81 int donor_idx = 0 ;
82 for (int ii=0; ii<d_len; ii++) {
83 if(isfinite(donor[ii])) {
84 donor_sites[donor_idx] = ii+1 ;
85 donor_idx++ ;
86 }
87 }
88
89 printf("end of donor/acceptor stuff\n");
90
91 Pre_score** matrices = new Pre_score*[nr_paths];
92 for (int z=0; z<nr_paths; z++)
93 matrices[z] = new Pre_score[dna_len * est_len];
94
95 char* est = (char*) malloc(sizeof(char)*est_len);
96 char* dna = (char*) malloc(sizeof(char)*dna_len);
97
98 for(idx=0;idx<est_len;idx++)
99 est[idx] = alphabet[lrand48()%4];
100
101 for(idx=0;idx<dna_len;idx++)
102 dna[idx] = alphabet[lrand48()%4];
103
104 double* prb = (double*) malloc(sizeof(double)*est_len);
105
106 for(idx=0;idx<est_len;idx++)
107 prb[idx] = lrand48() % 45;
108
109 double* matchmatrix = (double*) malloc(sizeof(double)*6);
110 matchmatrix[0] = 0;
111 matchmatrix[1] = -1;
112 matchmatrix[2] = -1;
113 matchmatrix[3] = -1;
114 matchmatrix[4] = -1;
115 matchmatrix[5] = 0;
116
117 REAL limits_array[10] = {19.999999999999993, 33.362010744001154, 55.651188044142465, 92.831776672255558, 154.85273653622525, 258.30993300297655, 430.88693800637651, 718.76273276092525, 1198.968500637882, 1999.9999999999982};
118 REAL pen_array[10] = {0.024414829630412388, 0.025836097651978609, 0.17473437525993854, 0.46902269254060858, -0.14383509580738071, -0.20766155250201096, -0.32734767088988909, -0.18507220858585502, -0.097862557903448902, -0.034886887426478781};
119
120 penalty_struct functions;
121
122 init_penalty_struct(functions);
123
124 functions.min_len = 20;
125 functions.max_len = 2000;
126 functions.len = 10;
127 functions.limits = &limits_array[0];
128 functions.penalties = &pen_array[0];
129
130 penalty_struct* qualityScores = (penalty_struct*) malloc(sizeof(penalty_struct)*30);
131
132 int num_params = 336;
133 double* params = (double*) malloc(sizeof(double)*num_params);
134
135 double parameter;
136 int ctr = 0;
137
138 string s;
139 ifstream infile("param.txt");
140
141 if (infile.fail())
142 printf("Could not open file param.txt for reading.\n");
143
144 while (getline(infile, s)) {
145 if (s.length() == 0 || s[0] == '#') continue;
146
147 istringstream iss(s);
148 iss >> parameter;
149 params[ctr] = parameter;
150 ctr++;
151 }
152
153 infile.close();
154
155 initialize_quality_plifs(qualityScores,params);
156
157 bool remove_duplicate_scores = true;
158 int* max_score_positions = new int[nr_paths*2];
159 mode currentMode = USE_QUALITY_SCORES;
160
161 printf("calling fill_matrix\n");
162 fill_matrix(nr_paths, matrices, est_len, dna_len, est, dna, prb, &functions, matchmatrix, qualityScores, donor, acceptor, remove_duplicate_scores, nr_donor_sites, donor_sites, max_score_positions);
163
164 printf("End of test_fm \n");
165 }