1 /* Last changed Time-stamp: <2008-10-09 16:23:36 ivo> */
19 static char rcsid[] = "$Id: findpath.c,v 1.2 2008/10/09 15:42:45 ivo Exp $";
22 #################################
24 #################################
28 #################################
30 #################################
32 PRIVATE const char *seq=NULL;
33 PRIVATE short *S=NULL, *S1=NULL;
35 PRIVATE move_t *path=NULL;
36 PRIVATE int path_fwd; /* 1: struc1->struc2, else struc2 -> struc1 */
40 /* NOTE: all variables are assumed to be uninitialized if they are declared as threadprivate
42 #pragma omp threadprivate(seq, S, S1, BP_dist, path, path_fwd)
47 #################################
48 # PRIVATE FUNCTION DECLARATIONS #
49 #################################
51 PRIVATE int *pair_table_to_loop_index (short *pt);
52 PRIVATE move_t *copy_moves(move_t *mvs);
53 PRIVATE int compare_ptable(const void *A, const void *B);
54 PRIVATE int compare_energy(const void *A, const void *B);
55 PRIVATE int compare_moves_when(const void *A, const void *B);
56 PRIVATE void free_intermediate(intermediate_t *i);
57 PRIVATE int find_path_once(const char *struc1, const char *struc2, int maxE, int maxl);
58 PRIVATE int try_moves(intermediate_t c, int maxE, intermediate_t *next, int dist);
61 #################################
62 # BEGIN OF FUNCTION DEFINITIONS #
63 #################################
65 PUBLIC void free_path(path_t *path){
68 while(tmp->s){ free(tmp->s); tmp++;}
73 PUBLIC int find_saddle(const char *sequence, const char *struc1, const char *struc2, int max) {
76 move_t *bestpath=NULL;
86 /* nummerically encode sequence */
87 S = encode_sequence(seq, 0);
88 S1 = encode_sequence(seq, 1);
94 if (maxl>max) maxl=max;
96 saddleE = find_path_once(struc1, struc2, maxE, maxl);
99 if (bestpath) free(bestpath);
104 free(path);path=NULL;
110 } while (maxl<2*max);
113 /* (re)set some globals */
119 PUBLIC void print_path(const char *seq, const char *struc) {
123 printf("%s\n%s %6.2f\n", seq, s, energy_of_structure(seq,s, 0));
124 qsort(path, BP_dist, sizeof(move_t), compare_moves_when);
125 for (d=0; d<BP_dist; d++) {
127 i = path[d].i; j=path[d].j;
128 if (i<0) { /* delete */
129 s[(-i)-1] = s[(-j)-1] = '.';
131 s[i-1] = '('; s[j-1] = ')';
133 printf("%s %6.2f - %6.2f\n", s, energy_of_structure(seq,s, 0), path[d].E/100.0);
138 PUBLIC path_t *get_path(const char *seq, const char *s1, const char* s2, int maxkeep) {
142 E = find_saddle(seq, s1, s2, maxkeep);
144 route = (path_t *)space((BP_dist+2)*sizeof(path_t));
146 qsort(path, BP_dist, sizeof(move_t), compare_moves_when);
149 /* memorize start of path */
150 route[0].s = strdup(s1);
151 route[0].en = energy_of_structure(seq, s1, 0);
153 for (d=0; d<BP_dist; d++) {
155 route[d+1].s = strdup(route[d].s);
156 i = path[d].i; j=path[d].j;
157 if (i<0) { /* delete */
158 route[d+1].s[(-i)-1] = route[d+1].s[(-j)-1] = '.';
160 route[d+1].s[i-1] = '('; route[d+1].s[j-1] = ')';
162 route[d+1].en = path[d].E/100.0;
166 /* memorize start of path */
168 route[BP_dist].s = strdup(s2);
169 route[BP_dist].en = energy_of_structure(seq, s2, 0);
171 for (d=0; d<BP_dist; d++) {
173 route[BP_dist-d-1].s = strdup(route[BP_dist-d].s);
176 if (i<0) { /* delete */
177 route[BP_dist-d-1].s[(-i)-1] = route[BP_dist-d-1].s[(-j)-1] = '.';
179 route[BP_dist-d-1].s[i-1] = '('; route[BP_dist-d-1].s[j-1] = ')';
181 route[BP_dist-d-1].en = path[d].E/100.0;
186 fprintf(stderr, "\n%s\n%s\n%s\n\n", seq, s1, s2);
187 for (d=0; d<=BP_dist; d++)
188 fprintf(stderr, "%s %6.2f\n", route[d].s, route[d].en);
189 fprintf(stderr, "%d\n", *num_entry);
192 free(path);path=NULL;
196 PRIVATE int try_moves(intermediate_t c, int maxE, intermediate_t *next, int dist) {
197 int *loopidx, len, num_next=0, en, oldE;
202 loopidx = pair_table_to_loop_index(c.pt);
204 for (mv=c.moves; mv->i!=0; mv++) {
206 if (mv->when>0) continue;
207 i = mv->i; j = mv->j;
208 pt = (short *) space(sizeof(short)*(len+1));
209 memcpy(pt, c.pt,(len+1)*sizeof(short));
210 if (j<0) { /*it's a delete move */
213 } else { /* insert move */
214 if ((loopidx[i] == loopidx[j]) && /* i and j belong to same loop */
215 (pt[i] == 0) && (pt[j]==0) /* ... and are unpaired */
221 continue; /* llegal move, try next; */
225 en = c.curr_en + energy_of_move_pt(c.pt, S, S1, i, j);
227 en = energy_of_structure_pt(seq, pt, S, S1, 0);
230 next[num_next].Sen = (en>oldE)?en:oldE;
231 next[num_next].curr_en = en;
232 next[num_next].pt = pt;
235 next[num_next++].moves = copy_moves(c.moves);
244 PRIVATE int find_path_once(const char *struc1, const char *struc2, int maxE, int maxl) {
247 int i, len, d, dist=0, result;
248 intermediate_t *current, *next;
250 pt1 = make_pair_table(struc1);
251 pt2 = make_pair_table(struc2);
252 len = (int) strlen(struc1);
254 mlist = (move_t *) space(sizeof(move_t)*len); /* bp_dist < n */
256 for (i=1; i<=len; i++) {
257 if (pt1[i] != pt2[i]) {
258 if (i<pt1[i]) { /* need to delete this pair */
260 mlist[dist].j = -pt1[i];
261 mlist[dist++].when = 0;
263 if (i<pt2[i]) { /* need to insert this pair */
265 mlist[dist].j = pt2[i];
266 mlist[dist++].when = 0;
272 current = (intermediate_t *) space(sizeof(intermediate_t)*(maxl+1));
274 current[0].Sen = current[0].curr_en = energy_of_structure_pt(seq, pt1, S, S1, 0);
275 current[0].moves = mlist;
276 next = (intermediate_t *) space(sizeof(intermediate_t)*(dist*maxl+1));
278 for (d=1; d<=dist; d++) { /* go through the distance classes */
279 int c, u, num_next=0;
282 for (c=0; current[c].pt != NULL; c++) {
283 num_next += try_moves(current[c], maxE, next+num_next, d);
286 for (cc=current; cc->pt != NULL; cc++) free_intermediate(cc);
287 current[0].Sen=INT_MAX;
290 /* remove duplicates via sort|uniq
291 if this becomes a bottleneck we can use a hash instead */
292 qsort(next, num_next, sizeof(intermediate_t),compare_ptable);
293 for (u=0,c=1; c<num_next; c++) {
294 if (memcmp(next[u].pt,next[c].pt,sizeof(short)*len)!=0) {
297 free_intermediate(next+c);
301 qsort(next, num_next, sizeof(intermediate_t),compare_energy);
302 /* free the old stuff */
303 for (cc=current; cc->pt != NULL; cc++) free_intermediate(cc);
304 for (u=0; u<maxl && u<num_next; u++) {
305 current[u] = next[u];
307 for (; u<num_next; u++)
308 free_intermediate(next+u);
312 path = current[0].moves;
313 result = current[0].Sen;
314 free(current[0].pt); free(current);
319 PRIVATE int *pair_table_to_loop_index (short *pt){
320 /* number each position by which loop it belongs to (positions start
328 stack = (int *) space(sizeof(int)*(length+1));
329 loop = (int *) space(sizeof(int)*(length+2));
332 for (i=1; i<=length; i++) {
333 if ((pt[i] != 0) && (i < pt[i])) { /* ( */
339 if ((pt[i] != 0) && (i > pt[i])) { /* ) */
342 l = loop[stack[hx-1]]; /* index of enclosing loop */
343 else l=0; /* external loop has index 0 */
345 nrerror("unbalanced brackets in make_pair_table");
352 #ifdef _DEBUG_LOOPIDX
353 fprintf(stderr,"begin loop index\n");
355 "....,....1....,....2....,....3....,....4"
356 "....,....5....,....6....,....7....,....8\n");
357 print_structure(pt, loop[0]);
358 for (i=1; i<=length; i++)
359 fprintf(stderr,"%2d ", loop[i]);
360 fprintf(stderr,"\n");
361 fprintf(stderr, "end loop index\n");
368 PRIVATE void free_intermediate(intermediate_t *i) {
376 PRIVATE int compare_ptable(const void *A, const void *B) {
377 intermediate_t *a, *b;
379 a = (intermediate_t *) A;
380 b = (intermediate_t *) B;
382 c = memcmp(a->pt, b->pt, a->pt[0]*sizeof(short));
384 if ((a->Sen - b->Sen) != 0) return (a->Sen - b->Sen);
385 return (a->curr_en - b->curr_en);
388 PRIVATE int compare_energy(const void *A, const void *B) {
389 intermediate_t *a, *b;
390 a = (intermediate_t *) A;
391 b = (intermediate_t *) B;
393 if ((a->Sen - b->Sen) != 0) return (a->Sen - b->Sen);
394 return (a->curr_en - b->curr_en);
397 PRIVATE int compare_moves_when(const void *A, const void *B) {
402 return(a->when - b->when);
405 PRIVATE move_t* copy_moves(move_t *mvs) {
407 new = (move_t *) space(sizeof(move_t)*(BP_dist+1));
408 memcpy(new,mvs,sizeof(move_t)*(BP_dist+1));
414 int main(int argc, char *argv[]) {
420 for (i=1; i<argc; i++) {
421 switch ( argv[i][1] ) {
422 case 'm': if (strcmp(argv[i],"-m")==0)
423 sscanf(argv[++i], "%d", &maxkeep);
425 case 'v': verbose = !strcmp(argv[i],"-v");
427 case 'd': if (strcmp(argv[i],"-d")==0)
428 sscanf(argv[++i], "%d", &dangles);
432 seq = get_line(stdin);
433 s1 = get_line(stdin);
434 s2 = get_line(stdin);
436 E = find_saddle(seq, s1, s2, maxkeep);
437 printf("saddle_energy = %6.2f\n", E/100.);
445 route = get_path(seq, s1, s2, maxkeep);
446 for (r=route; r->s; r++) {
447 printf("%s %6.2f - %6.2f\n", r->s, energy_of_struct(seq,r->s), r->en);
451 free(seq); free(s1); free(s2);
452 return(EXIT_SUCCESS);