1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include <limits.h>
6 #include <time.h>
7
8 #include "move_set_pk.h"
9
10 extern "C" {
11 #include "pair_mat.h"
12 }
13
14 /* maximum degeneracy value - if degeneracy is greater than this, program segfaults*/
15 #define MAX_DEGEN 100
16 #define MINGAP 3
17
18 //#define space(a) malloc(a)
19
20 #define bool int
21 #define true 1
22 #define false 0
23
compare(Structure * lhs,Structure * rhs)24 int compare(Structure *lhs, Structure *rhs)
25 {
26 return (*lhs)<(*rhs);
27 }
28
find_min(Structure * arr[MAX_DEGEN],int begin,int end)29 int find_min(Structure *arr[MAX_DEGEN], int begin, int end) {
30 Structure *min = arr[begin];
31 int min_num = begin;
32 int i;
33
34 for (i=begin+1; i<end; i++) {
35 if (compare(arr[i], min)) {
36 min = arr[i];
37 min_num = i;
38 }
39 }
40 return min_num;
41 }
42
equals(const Structure * first,const Structure * second)43 int equals(const Structure *first, const Structure *second)
44 {
45 return (*first)==(*second);
46 }
47
48 /* ############################## DECLARATION #####################################*/
49 /* private functions & declarations*/
50
51 static int cnt_move = 0;
count_move()52 int count_move() {return cnt_move;}
53
54 void print_str_pk(FILE *out, short *str);
55
56 /*check if base is lone*/
57 int lone_base(short *pt, int i);
58
59 /* internal struct with moves, sequence, degeneracy and options*/
60 typedef struct _Encoded {
61 // sequence
62 const char *seq;
63 short *s0;
64 short *s1;
65
66 /* moves*/
67 int bp_left;
68 int bp_right;
69
70 /* options*/
71 //int noLP;
72 int verbose_lvl;
73 int first;
74 int shift;
75 int all_neighs; // sould be on for shifts!
76
77 /* degeneracy*/
78 int begin_unpr;
79 int begin_pr;
80 int end_unpr;
81 int end_pr;
82 Structure *processed[MAX_DEGEN];
83 Structure *unprocessed[MAX_DEGEN];
84 int current_en;
85
86 /* moves in random (needs to be freed afterwards)*/
87 int *moves_from;
88 int *moves_to;
89 int num_moves;
90
91 /* function for flooding */
92 int (*funct) (Structure*, Structure*);
93 } Encoded;
94
95 /* frees all things allocated by degeneracy...*/
free_degen(Encoded * Enc)96 void free_degen(Encoded *Enc)
97 {
98 int i;
99 for (i=Enc->begin_unpr; i<Enc->end_unpr; i++) {
100 if (Enc->unprocessed[i]) {
101 delete Enc->unprocessed[i];
102 Enc->unprocessed[i]=NULL;
103 }
104 }
105 for (i=Enc->begin_pr; i<Enc->end_pr; i++) {
106 if (Enc->processed[i]) {
107 delete Enc->processed[i];
108 Enc->processed[i]=NULL;
109 }
110 }
111 Enc->begin_pr=0;
112 Enc->begin_unpr=0;
113 Enc->end_pr=0;
114 Enc->end_unpr=0;
115 }
116
117 /* ############################## IMPLEMENTATION #####################################*/
118
119
print_str_pk(FILE * out,short * str)120 void print_str_pk(FILE *out, short *str)
121 {
122 fprintf(out, "%s", pt_to_str_pk(str).c_str());
123 }
124
do_move(short * pt,int bp_left,int bp_right)125 inline void do_move(short *pt, int bp_left, int bp_right)
126 {
127 /* delete*/
128 if (bp_left<0) {
129 pt[-bp_left]=0;
130 pt[-bp_right]=0;
131 } else { /* insert*/
132 pt[bp_left]=bp_right;
133 pt[bp_right]=bp_left;
134 }
135 }
136
137 /* done with all structures along the way to deepest*/
update_deepest(Encoded * Enc,Structure * str,Structure * min)138 int update_deepest(Encoded *Enc, Structure *str, Structure *min)
139 {
140 /* apply move + get its energy*/
141 int tmp_en;
142 tmp_en = str->energy + str->MakeMove(Enc->seq, Enc->s0, Enc->s1, Enc->bp_left, Enc->bp_right);
143
144 /* use f_point if we have it */
145 if (Enc->funct) {
146 int end = Enc->funct(str, min);
147
148 // undo moves
149 str->UndoMove();
150 Enc->bp_left=0;
151 Enc->bp_right=0;
152
153 return (end?1:0);
154 }
155
156 if (Enc->verbose_lvl>1) { fprintf(stderr, " "); print_str_pk(stderr, str->str); fprintf(stderr, " %d\n", tmp_en); }
157
158 /* better deepest*/
159 if (tmp_en < min->energy) {
160
161 *min = *str;
162
163 /* delete degeneracy*/
164 free_degen(Enc);
165
166 /* undo moves*/
167 str->UndoMove();
168 Enc->bp_left=0;
169 Enc->bp_right=0;
170 return 1;
171 }
172
173 /* degeneracy*/
174 if ((str->energy == min->energy) && (Enc->current_en == min->energy)) {
175 int found = 0;
176 int i;
177 for (i=Enc->begin_pr; i<Enc->end_pr; i++) {
178 if (equals(Enc->processed[i], str)) {
179 found = 1;
180 break;
181 }
182 }
183 for (i=Enc->begin_unpr; !found && i<Enc->end_unpr; i++) {
184 if (equals(Enc->unprocessed[i], str)) {
185 found = 1;
186 break;
187 }
188 }
189
190 if (!found) {
191 //fprintf(stderr, "%s %6.2f\n", pt_to_str_pk(str->str).c_str(), str->energy);
192 Enc->unprocessed[Enc->end_unpr] = new Structure(*str);
193 if (Enc->end_unpr == MAX_DEGEN) {
194 fprintf(stderr, "ERROR: Degeneracy too high!!! %s\n", Enc->seq);
195 for (int j=Enc->begin_unpr; j<Enc->end_unpr; j++) {
196 fprintf(stderr, "D%3d %s %6.2f\n", j, pt_to_str_pk(Enc->unprocessed[j]->str).c_str(), Enc->unprocessed[j]->energy/100.0);
197 }
198 fprintf(stderr, "\n");
199 for (int j=Enc->begin_pr; j<Enc->end_pr; j++) {
200 fprintf(stderr, "D%3d %s %6.2f\n", j, pt_to_str_pk(Enc->processed[j]->str).c_str(), Enc->processed[j]->energy/100.0);
201 }
202 exit(EXIT_FAILURE);
203 }
204 Enc->end_unpr++;
205 }
206 }
207
208 /* undo moves*/
209 str->UndoMove();
210 Enc->bp_left=0;
211 Enc->bp_right=0;
212 return 0;
213 }
214
215
216 /* deletions move set*/
deletions_pk(Encoded * Enc,Structure * str,Structure * minim)217 int deletions_pk(Encoded *Enc, Structure *str, Structure *minim)
218 {
219 int cnt = 0;
220 short *pt = str->str;
221 int len = pt[0];
222 int i;
223
224 for (i=1; i<=len; i++) {
225 if (pt[i] && pt[i]>i) { /* '('*/
226 Enc->bp_left=-i;
227 Enc->bp_right=-pt[i];
228
229
230 cnt += update_deepest(Enc, str, minim);
231 /* in case useFirst is on and structure is found, end*/
232 if (Enc->first && cnt > 0) return cnt;
233 }
234 }
235 return cnt;
236 }
237
238 /* compatible base pair?*/
compat(char a,char b)239 inline bool compat(char a, char b) {
240 if (a=='A' && b=='U') return true;
241 if (a=='C' && b=='G') return true;
242 if (a=='G' && b=='U') return true;
243 if (a=='U' && b=='A') return true;
244 if (a=='G' && b=='C') return true;
245 if (a=='U' && b=='G') return true;
246 /* and with T's*/
247 if (a=='A' && b=='T') return true;
248 if (a=='T' && b=='A') return true;
249 if (a=='G' && b=='T') return true;
250 if (a=='T' && b=='G') return true;
251 return false;
252 }
253
254 /* try insert base pair (i,j)*/
try_insert(const short * pt,const char * seq,int i,int j)255 inline bool try_insert(const short *pt, const char *seq, int i, int j)
256 {
257 if (i<=0 || j<=0 || i>pt[0] || j>pt[0]) return false;
258 return (j-i>MINGAP && pt[j]==0 && pt[i]==0 && compat(seq[i-1], seq[j-1]));
259 }
260
261 /* try switch base pair (i,j)*/
try_shift(const short * pt,const char * seq,int i,int j)262 inline bool try_shift(const short *pt, const char *seq, int i, int j)
263 {
264 if (i<=0 || j<=0 || i>pt[0] || j>pt[0]) return false;
265 return (j-i>MINGAP && compat(seq[i-1], seq[j-1]) &&
266 ((pt[i]==0 && pt[j]!=0) || (pt[i]!=0 && pt[j]==0)) );
267 }
268
269 // try insert base pair (i,j)
try_insert_seq(const char * seq,int i,int j)270 inline bool try_insert_seq(const char *seq, int i, int j)
271 {
272 if (i<=0 || j<=0) return false;
273 return (j-i>MINGAP && compat(seq[i-1], seq[j-1]));
274 }
275
276 /* insertions move set*/
insertions_pk(Encoded * Enc,Structure * str,Structure * minim)277 int insertions_pk(Encoded *Enc, Structure *str, Structure *minim)
278 {
279 int cnt = 0;
280 short *pt = str->str;
281 int len = pt[0];
282 int i,j;
283
284 for (i=1; i<=len; i++) {
285 if (pt[i]==0) {
286 for (j=i+1; j<=len; j++) {
287 /* end if found closing bracket*/
288 //if (pt[j]!=0 && pt[j]<j) break; //')'
289 /*if (pt[j]!=0 && pt[j]>j) { //'('
290 j = pt[j];
291 continue;
292 }*/
293 /* if conditions are met, do insert*/
294 if (try_insert(pt, Enc->seq, i, j)) {
295 INS_FLAG iflag = str->CanInsert(i, j);
296 if (iflag != NO_INS && (Enc->all_neighs || iflag == REG_INS || iflag == INSIDE_PK)) {
297 Enc->bp_left=i;
298 Enc->bp_right=j;
299
300 //fprintf(stderr, "CI: %4d %4d\n", i, j);
301
302
303 cnt += update_deepest(Enc, str, minim);
304 /* in case useFirst is on and structure is found, end*/
305 if (Enc->first && cnt > 0) return cnt;
306 } /*else {
307 str->str[i] = j;
308 str->str[j] = i;
309 fprintf(stdout, "%s BAD\n", pt_to_str_pk(str->str).c_str());
310 str->str[i] = 0;
311 str->str[j] = 0;
312 }*/
313 }
314 }
315 }
316 }
317 return cnt;
318 }
319
320 /* shifts move set*/
shifts_pk(Encoded * Enc,Structure * str,Structure * minim)321 int shifts_pk(Encoded *Enc, Structure *str, Structure *minim)
322 {
323 int cnt = 0;
324 short *pt = str->str;
325 int len = pt[0];
326 int i,j;
327
328 for (i=1; i<=len; i++) {
329 if (pt[i]==0) continue;
330
331 // first '('
332 if (pt[i]>i) {
333
334 for (j=i+1; j<=len; j++) {
335
336 if (try_shift(pt, Enc->seq, i, j)) {
337 INS_FLAG iflag = str->CanShift(i, j);
338 if (iflag != NO_INS && (Enc->all_neighs || iflag == REG_INS || iflag == INSIDE_PK)) {
339 Enc->bp_left=i;
340 Enc->bp_right=j;
341
342 //fprintf(stderr, "CS: %4d %4d\n", i, j);
343
344
345 cnt += update_deepest(Enc, str, minim);
346 /* in case useFirst is on and structure is found, end*/
347 if (Enc->first && cnt > 0) return cnt;
348 }
349 }
350 }
351 }
352
353 // then ')'
354 if (pt[i]<i) {
355 for (j=1; j<i; j++) {
356
357 if (try_shift(pt, Enc->seq, i, j)) {
358 INS_FLAG iflag = str->CanShift(i, j);
359 if (iflag != NO_INS && (Enc->all_neighs || iflag == REG_INS || iflag == INSIDE_PK)) {
360 Enc->bp_left=i;
361 Enc->bp_right=j;
362
363 //fprintf(stderr, "CS: %4d %4d\n", i, j);
364
365
366 cnt += update_deepest(Enc, str, minim);
367 /* in case useFirst is on and structure is found, end*/
368 if (Enc->first && cnt > 0) return cnt;
369 }
370 }
371 }
372 }
373
374 }
375 return cnt;
376 }
377
378 /* move to deepest (or first) neighbour*/
move_set(Encoded * Enc,Structure * str_in)379 int move_set(Encoded *Enc, Structure *str_in)
380 {
381 /* count how many times called*/
382 cnt_move++;
383
384 /* count better neighbours*/
385 int cnt = 0;
386
387 /* deepest descent*/
388 Structure *str = new Structure(*str_in);
389 Structure *min = new Structure(*str);
390 Enc->current_en = str->energy;
391
392 if (Enc->verbose_lvl>1) { fprintf(stderr, " start of MS:\n "); print_str_pk(stderr, str->str); fprintf(stderr, " %d\n\n", str->energy); }
393
394 /* if using first dont do all of them*/
395 bool end = false;
396 /* insertions*/
397 if (!end) cnt += insertions_pk(Enc, str, min);
398 if (Enc->first && cnt>0) end = true;
399 if (Enc->verbose_lvl>1) fprintf(stderr, "\n");
400
401 /* deletions*/
402 if (!end) cnt += deletions_pk(Enc, str, min);
403 if (Enc->first && cnt>0) end = true;
404
405 /* shifts*/
406 if (Enc->shift) {
407 if (!end) cnt += shifts_pk(Enc, str, min);
408 if (Enc->first && cnt>0) end = true;
409 }
410
411 /* if degeneracy occurs, solve it!*/
412 if (!end && (Enc->end_unpr - Enc->begin_unpr)>0) {
413 Enc->processed[Enc->end_pr] = str;
414 Enc->end_pr++;
415 str = Enc->unprocessed[Enc->begin_unpr];
416 Enc->unprocessed[Enc->begin_unpr]=NULL;
417 Enc->begin_unpr++;
418 delete min;
419 cnt += move_set(Enc, str);
420 } else {
421 /* write output to str*/
422 *str = *min;
423 delete min;
424 }
425
426 /* resolve degeneracy in local minima*/
427 if ((Enc->end_pr - Enc->begin_pr)>0) {
428 Enc->processed[Enc->end_pr]=str;
429 Enc->end_pr++;
430
431 int min = find_min(Enc->processed, Enc->begin_pr, Enc->end_pr);
432 Structure *tmp = Enc->processed[min];
433 Enc->processed[min] = Enc->processed[Enc->begin_pr];
434 Enc->processed[Enc->begin_pr] = tmp;
435 str = Enc->processed[Enc->begin_pr];
436 Enc->begin_pr++;
437 free_degen(Enc);
438 }
439
440 if (Enc->verbose_lvl>1 && !(Enc->first)) { fprintf(stderr, "\n end of MS:\n "); print_str_pk(stderr, str->str); fprintf(stderr, " %d\n\n", str->energy); }
441
442 *str_in = *str;
443 delete str;
444
445 return cnt;
446 }
447
construct_moves(Encoded * Enc,short * structure)448 void construct_moves(Encoded *Enc, short *structure)
449 {
450 /* generate all possible moves (less than n^2)*/
451 Enc->num_moves = 0;
452 int i;
453 for (i=1; i<=structure[0]; i++) {
454 if (structure[i]!=0) {
455 if (structure[i]<i) continue;
456 Enc->moves_from[Enc->num_moves]=-i;
457 Enc->moves_to[Enc->num_moves]=-structure[i];
458 Enc->num_moves++;
459 //fprintf(stderr, "add d(%d, %d)\n", i, str.structure[i]);
460 } else {
461 int j;
462 for (j=i+1; j<=structure[0]; j++) {
463 //fprintf(stderr, "check (%d, %d)\n", i, j);
464 if (structure[j]==0) {
465 if (try_insert_seq(Enc->seq,i,j)) {
466 Enc->moves_from[Enc->num_moves]=i;
467 Enc->moves_to[Enc->num_moves]=j;
468 Enc->num_moves++;
469 //fprintf(stderr, "add i(%d, %d)\n", i, j);
470 continue;
471 }
472 } /*else if (structure[j]>j) { // '('
473 j = structure[j];
474 } else break;*/
475 }
476 }
477 }
478
479 /* permute them */
480 for (i=0; i<Enc->num_moves-1; i++) {
481 int rnd = rand();
482 rnd = rnd % (Enc->num_moves-i) + i;
483 int swp;
484 swp = Enc->moves_from[i];
485 Enc->moves_from[i]=Enc->moves_from[rnd];
486 Enc->moves_from[rnd]=swp;
487 swp = Enc->moves_to[i];
488 Enc->moves_to[i]=Enc->moves_to[rnd];
489 Enc->moves_to[rnd]=swp;
490 }
491 }
492
move_rset(Encoded * Enc,Structure * str_in)493 int move_rset(Encoded *Enc, Structure *str_in)
494 {
495 /* count how many times called*/
496 cnt_move++;
497
498 /* count better neighbours*/
499 int cnt = 0;
500
501 /* deepest descent*/
502 Structure *str = new Structure(*str_in);
503 Structure *min = new Structure(*str);
504 Enc->current_en = str->energy;
505
506 if (Enc->verbose_lvl>1) { fprintf(stderr, " start of MR:\n "); print_str_pk(stderr, str->str); fprintf(stderr, " %d\n\n", str->energy); }
507
508 // construct and permute possible moves
509 construct_moves(Enc, str->str);
510
511 /* find first the lower one*/
512 int i;
513 for (i=0; i<Enc->num_moves; i++) {
514 Enc->bp_left = Enc->moves_from[i];
515 Enc->bp_right = Enc->moves_to[i];
516 cnt = update_deepest(Enc, str, min);
517 if (cnt) break;
518 }
519
520 /* if degeneracy occurs, solve it!*/
521 if (!cnt && (Enc->end_unpr - Enc->begin_unpr)>0) {
522 Enc->processed[Enc->end_pr] = str;
523 Enc->end_pr++;
524 str = Enc->unprocessed[Enc->begin_unpr];
525 Enc->unprocessed[Enc->begin_unpr]=NULL;
526 Enc->begin_unpr++;
527 delete min;
528 cnt += move_rset(Enc, str);
529 } else {
530 /* write output to str*/
531
532 *str = *min;
533 delete min;
534 }
535
536 /* resolve degeneracy in local minima*/
537 if ((Enc->end_pr - Enc->begin_pr)>0) {
538 Enc->processed[Enc->end_pr]=str;
539 Enc->end_pr++;
540
541 int min = find_min(Enc->processed, Enc->begin_pr, Enc->end_pr);
542 Structure *tmp = Enc->processed[min];
543 Enc->processed[min] = Enc->processed[Enc->begin_pr];
544 Enc->processed[Enc->begin_pr] = tmp;
545 str = Enc->processed[Enc->begin_pr];
546 Enc->begin_pr++;
547 free_degen(Enc);
548 }
549
550 *str_in = *str;
551 delete str;
552
553 return cnt;
554 }
555
556 /*check if base is lone*/
557 /*int lone_base(short *pt, int i)
558 {
559 if (i<=0 || i>pt[0]) return 0;
560 // is not a base pair
561 if (pt[i]==0) return 0;
562
563 // base is lone:
564 if (i-1>0) {
565 // is base pair and is the same bracket
566 if (pt[i-1]!=0 && ((pt[i-1]<pt[pt[i-1]]) == (pt[i]<pt[pt[i]]))) return 0;
567 }
568
569 if (i+1<=pt[0]) {
570 if (pt[i+1]!=0 && ((pt[i-1]<pt[pt[i-1]]) == (pt[i]<pt[pt[i]]))) return 0;
571 }
572
573 return 1;
574 }*/
575
move_standard_pk_pt(const char * seq,Structure * str,short * s0,short * s1,enum MOVE_TYPE type,int shifts,int verbosity_level)576 int move_standard_pk_pt(const char *seq,
577 Structure *str,
578 short *s0,
579 short *s1,
580 enum MOVE_TYPE type,
581 int shifts,
582 int verbosity_level)
583 {
584 switch (type){
585 case GRADIENT: move_gradient_pk(seq, str, s0, s1, shifts, verbosity_level); break;
586 case FIRST: move_first_pk(seq, str, s0, s1, shifts, verbosity_level); break;
587 case ADAPTIVE: move_adaptive_pk(seq, str, s0, s1, shifts, verbosity_level); break;
588 }
589
590 return str->energy;
591 }
592
move_standard_pk(const char * seq,char * struc,enum MOVE_TYPE type,int shifts,int verbosity_level)593 int move_standard_pk(const char *seq,
594 char *struc,
595 enum MOVE_TYPE type,
596 int shifts,
597 int verbosity_level)
598 {
599 make_pair_matrix();
600 short *s0 = encode_sequence(seq, 0);
601 short *s1 = encode_sequence(seq, 1);
602
603 Structure *str = new Structure(seq, struc, s0, s1);
604
605 int energy = move_standard_pk_pt(seq, str, s0, s1, type, shifts, verbosity_level);
606
607 free(s0);
608 free(s1);
609
610 pt_to_chars_pk(str->str, struc);
611 delete str;
612
613 return energy;
614 }
615
move_gradient_pk(const char * seq,Structure * str,short * s0,short * s1,int shifts,int verbosity_level)616 int move_gradient_pk(const char *seq,
617 Structure *str,
618 short *s0,
619 short *s1,
620 int shifts,
621 int verbosity_level)
622 {
623 cnt_move = 0;
624
625 Encoded enc;
626 enc.seq = seq;
627 enc.s0 = s0;
628 enc.s1 = s1;
629
630 /* moves*/
631 enc.bp_left=0;
632 enc.bp_right=0;
633
634 /* options*/
635 enc.verbose_lvl=verbosity_level;
636 enc.first=0;
637 enc.all_neighs=shifts;
638 enc.shift = shifts;
639
640 /* degeneracy*/
641 enc.begin_unpr=0;
642 enc.begin_pr=0;
643 enc.end_unpr=0;
644 enc.end_pr=0;
645 enc.current_en=0;
646
647 // function
648 enc.funct=NULL;
649
650 int i;
651 for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
652
653 /*struct_en str;
654 str.structure = allocopy(ptable);
655 str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);
656 */
657
658 while (move_set(&enc, str)!=0) {
659 free_degen(&enc);
660 }
661 free_degen(&enc);
662
663 return str->energy;
664 }
665
move_first_pk(const char * seq,Structure * str,short * s0,short * s1,int shifts,int verbosity_level)666 int move_first_pk(const char *seq,
667 Structure *str,
668 short *s0,
669 short *s1,
670 int shifts,
671 int verbosity_level)
672 {
673 cnt_move = 0;
674
675 Encoded enc;
676 enc.seq = seq;
677 enc.s0 = s0;
678 enc.s1 = s1;
679
680 /* moves*/
681 enc.bp_left=0;
682 enc.bp_right=0;
683
684 /* options*/
685 enc.verbose_lvl=verbosity_level;
686 enc.first=1;
687 enc.all_neighs=shifts;
688 enc.shift = shifts;
689
690 /* degeneracy*/
691 enc.begin_unpr=0;
692 enc.begin_pr=0;
693 enc.end_unpr=0;
694 enc.end_pr=0;
695 enc.current_en=0;
696
697 // function
698 enc.funct=NULL;
699
700 int i;
701 for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
702
703 while (move_set(&enc, str)!=0) {
704 free_degen(&enc);
705 }
706 free_degen(&enc);
707
708 return str->energy;
709 }
710
move_adaptive_pk(const char * seq,Structure * str,short * s0,short * s1,int shifts,int verbosity_level)711 int move_adaptive_pk(const char *seq,
712 Structure *str,
713 short *s0,
714 short *s1,
715 int shifts,
716 int verbosity_level)
717 {
718 srand(time(NULL));
719
720 cnt_move = 0;
721
722 Encoded enc;
723 enc.seq = seq;
724 enc.s0 = s0;
725 enc.s1 = s1;
726
727 /* moves*/
728 enc.bp_left=0;
729 enc.bp_right=0;
730
731 /* options*/
732 enc.verbose_lvl=verbosity_level;
733 enc.first=1;
734 enc.all_neighs=shifts;
735 enc.shift = shifts;
736
737 /* degeneracy*/
738 enc.begin_unpr=0;
739 enc.begin_pr=0;
740 enc.end_unpr=0;
741 enc.end_pr=0;
742 enc.current_en=0;
743
744 // function
745 enc.funct=NULL;
746
747 // allocate memory for moves
748 int length = str->str[0];
749 enc.moves_from = (int*) space(length*length*sizeof(int));
750 enc.moves_to = (int*) space(length*length*sizeof(int));
751
752 int i;
753 for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
754
755 while (move_rset(&enc, str)!=0) {
756 free_degen(&enc);
757 }
758 free_degen(&enc);
759
760 free(enc.moves_from);
761 free(enc.moves_to);
762
763 return str->energy;
764 }
765
browse_neighs_pk(const char * seq,char * struc,int shifts,int verbosity_level,int (* funct)(Structure *,Structure *))766 int browse_neighs_pk( const char *seq,
767 char *struc,
768 int shifts,
769 int verbosity_level,
770 int (*funct) (Structure*, Structure*))
771
772 {
773 make_pair_matrix();
774 short *s0 = encode_sequence(seq, 0);
775 short *s1 = encode_sequence(seq, 1);
776
777 Structure *str = new Structure(seq, struc, s0, s1);
778
779 int res = browse_neighs_pk_pt(seq, str, s0, s1, shifts, verbosity_level, funct);
780
781 free(s0);
782 free(s1);
783
784 delete str;
785
786 return res;
787 }
788
browse_neighs_pk_pt(const char * seq,Structure * str,short * s0,short * s1,int shifts,int verbosity_level,int (* funct)(Structure *,Structure *))789 int browse_neighs_pk_pt(const char *seq,
790 Structure *str,
791 short *s0,
792 short *s1,
793 int shifts,
794 int verbosity_level,
795 int (*funct) (Structure*, Structure*))
796 {
797 cnt_move = 0;
798
799 Encoded enc;
800 enc.seq = seq;
801 enc.s0 = s0;
802 enc.s1 = s1;
803
804 /* moves*/
805 enc.bp_left=0;
806 enc.bp_right=0;
807
808 /* options*/
809 enc.verbose_lvl=verbosity_level;
810 enc.first=1;
811 enc.all_neighs=1;
812 enc.shift = shifts;
813
814 /* degeneracy*/
815 enc.begin_unpr=0;
816 enc.begin_pr=0;
817 enc.end_unpr=0;
818 enc.end_pr=0;
819 enc.current_en=0;
820
821 // function
822 enc.funct=funct;
823
824 int i;
825 for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
826
827 move_set(&enc, str);
828 free_degen(&enc);
829
830 return str->energy;
831 }
832
833 /*
834 // sample usage:
835 int main() {
836 char seq[20] = "ACCCCCCTCTGTAGGGGGA";
837 char str[20] = ".((.(.........).)).";
838
839 // move to the local minimum and display it
840 int energy = move_standard(seq, str, GRADIENT, 0, 0, 0);
841 fprintf(stdout, "%s %6.2f\n\n", str, energy/100.0);
842
843 //now create an array of every structure in neighbourhood of str structure
844 struct_en *list = NULL;
845 int list_length = 0;
846
847 int get_list(struct_en *new_one, struct_en *old_one)
848 {
849 // enlarge the list
850 list_length++;
851 list = (struct_en*) realloc(list, list_length*sizeof(struct_en));
852
853 // copy the structure
854 list[list_length-1].energy = new_one->energy;
855 list[list_length-1].structure = allocopy(new_one->structure);
856
857 // we want to continue -> return 0
858 return 0;
859 }
860 browse_neighs(seq, str, 0, 0, 0, get_list);
861
862 // print them and free the memory:
863 int i;
864 for (i=0; i<list_length; i++) {
865 print_stren(stdout, &list[i]);
866 free(list[i].structure);
867 }
868 free(list);
869
870 return 0;
871 }*/
872
873
874