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