/* gcc -fopenmp -g3 -DTEST_FINDPATH findpath.c -o FINDpath -lRNA -lm -I ../H/ -L ./ */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include "ViennaRNA/findpath.h"
#include "ViennaRNA/data_structures.h"
#include "ViennaRNA/model.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/fold.h"
#include "ViennaRNA/cofold.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/utils.h"
#include "ViennaRNA/structure_utils.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#define LOOP_EN
/**
* @brief
*/
typedef struct move {
int i; /* i,j>0 insert; i,j<0 delete */
int j;
int when; /* 0 if still available, else resulting distance from start */
int E;
} move_t;
/**
* @brief
*/
typedef struct intermediate {
short *pt; /**< @brief pair table */
int Sen; /**< @brief saddle energy so far */
int curr_en; /**< @brief current energy */
move_t *moves; /**< @brief remaining moves to target */
} intermediate_t;
/*
#################################
# GLOBAL VARIABLES #
#################################
*/
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
PRIVATE int BP_dist;
PRIVATE move_t *path=NULL;
PRIVATE int path_fwd; /* 1: struc1->struc2, else struc2 -> struc1 */
PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
#ifdef _OPENMP
/* NOTE: all variables are assumed to be uninitialized if they are declared as threadprivate
*/
#pragma omp threadprivate(BP_dist, path, path_fwd, backward_compat_compound)
#endif
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE move_t *copy_moves(move_t *mvs);
PRIVATE int compare_ptable(const void *A, const void *B);
PRIVATE int compare_energy(const void *A, const void *B);
PRIVATE int compare_moves_when(const void *A, const void *B);
PRIVATE void free_intermediate(intermediate_t *i);
#ifdef TEST_FINDPATH
/* TEST_FINDPATH, COFOLD */
PRIVATE void usage(void);
#endif
PRIVATE int find_path_once(vrna_fold_compound_t *vc, const char *struc1, const char *struc2, int maxE, int maxl);
PRIVATE int try_moves(vrna_fold_compound_t *vc, intermediate_t c, int maxE, intermediate_t *next, int dist);
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
PUBLIC void free_path(vrna_path_t *path){
vrna_path_t *tmp = path;
if(tmp){
while(tmp->s){ free(tmp->s); tmp++;}
free(path);
}
}
PUBLIC int
find_saddle(const char *seq,
const char *struc1,
const char *struc2,
int max){
int maxE;
char *sequence;
vrna_fold_compound_t *vc;
vrna_md_t md;
vc = NULL;
set_model_details(&md);
if(backward_compat_compound){
if(!strcmp(seq, backward_compat_compound->sequence)){ /* check if sequence is the same as before */
md.window_size = backward_compat_compound->length;
if(!memcmp(&md, &(backward_compat_compound->params->model_details), sizeof(vrna_md_t))){ /* check if model_details are the same as before */
vc = backward_compat_compound; /* re-use previous vrna_fold_compound_t */
}
}
}
if(!vc){
vrna_fold_compound_free(backward_compat_compound);
sequence = vrna_cut_point_insert(seq, cut_point);
backward_compat_compound = vc = vrna_fold_compound(sequence, &md, VRNA_OPTION_EVAL_ONLY);
free(sequence);
}
maxE = vrna_path_findpath_saddle(vc, struc1, struc2, max);
return maxE;
}
PUBLIC int
vrna_path_findpath_saddle( vrna_fold_compound_t *vc,
const char *struc1,
const char *struc2,
int max) {
int maxl, maxE;
const char *tmp;
move_t *bestpath=NULL;
int dir;
path_fwd = 0;
maxE = INT_MAX - 1;
maxl = 1;
do {
int saddleE;
path_fwd = !path_fwd;
if (maxl>max) maxl=max;
if(path) free(path);
saddleE = find_path_once(vc, struc1, struc2, maxE, maxl);
if (saddleE<maxE) {
maxE = saddleE;
if (bestpath) free(bestpath);
bestpath = path;
path = NULL;
dir = path_fwd;
} else{
free(path);path=NULL;
}
tmp=struc1;
struc1=struc2;
struc2=tmp;
maxl *=2;
} while (maxl<2*max);
/* (re)set some globals */
path=bestpath;
path_fwd = dir;
return maxE;
}
PUBLIC vrna_path_t *
get_path( const char *seq,
const char *s1,
const char* s2,
int maxkeep){
vrna_path_t *route = NULL;
char *sequence = NULL;
vrna_fold_compound_t *vc = NULL;
vrna_md_t md;
set_model_details(&md);
if(backward_compat_compound){
if(!strcmp(seq, backward_compat_compound->sequence)){ /* check if sequence is the same as before */
if(!memcmp(&md, &(backward_compat_compound->params->model_details), sizeof(vrna_md_t))){ /* check if model_details are the same as before */
vc = backward_compat_compound; /* re-use previous vrna_fold_compound_t */
}
}
}
if(!vc){
vrna_fold_compound_free(backward_compat_compound);
sequence = vrna_cut_point_insert(seq, cut_point);
backward_compat_compound = vc = vrna_fold_compound(sequence, &md, VRNA_OPTION_EVAL_ONLY);
free(sequence);
}
route = vrna_path_findpath(vc, s1, s2, maxkeep);
return (route);
}
PUBLIC vrna_path_t *
vrna_path_findpath(vrna_fold_compound_t *vc,
const char *s1,
const char* s2,
int maxkeep){
int E, d;
vrna_path_t *route=NULL;
E = vrna_path_findpath_saddle(vc, s1, s2, maxkeep);
route = (vrna_path_t *)vrna_alloc((BP_dist+2)*sizeof(vrna_path_t));
qsort(path, BP_dist, sizeof(move_t), compare_moves_when);
if (path_fwd) {
/* memorize start of path */
route[0].s = strdup(s1);
route[0].en = vrna_eval_structure(vc, s1);
for (d=0; d<BP_dist; d++) {
int i,j;
route[d+1].s = strdup(route[d].s);
i = path[d].i; j=path[d].j;
if (i<0) { /* delete */
route[d+1].s[(-i)-1] = route[d+1].s[(-j)-1] = '.';
} else {
route[d+1].s[i-1] = '('; route[d+1].s[j-1] = ')';
}
route[d+1].en = path[d].E/100.0;
}
}
else {
/* memorize start of path */
route[BP_dist].s = strdup(s2);
route[BP_dist].en = vrna_eval_structure(vc, s2);
for (d=0; d<BP_dist; d++) {
int i,j;
route[BP_dist-d-1].s = strdup(route[BP_dist-d].s);
i = path[d].i;
j = path[d].j;
if (i<0) { /* delete */
route[BP_dist-d-1].s[(-i)-1] = route[BP_dist-d-1].s[(-j)-1] = '.';
} else {
route[BP_dist-d-1].s[i-1] = '('; route[BP_dist-d-1].s[j-1] = ')';
}
route[BP_dist-d-1].en = path[d].E/100.0;
}
}
#if _DEBUG_FINDPATH_
fprintf(stderr, "\n%s\n%s\n%s\n\n", seq, s1, s2);
for (d=0; d<=BP_dist; d++)
fprintf(stderr, "%s %6.2f\n", route[d].s, route[d].en);
fprintf(stderr, "%d\n", *num_entry);
#endif
free(path);path=NULL;
return (route);
}
PRIVATE int
try_moves(vrna_fold_compound_t *vc,
intermediate_t c,
int maxE,
intermediate_t *next,
int dist){
int *loopidx, len, num_next=0, en, oldE;
move_t *mv;
short *pt;
len = c.pt[0];
loopidx = vrna_loopidx_from_ptable(c.pt);
oldE = c.Sen;
for (mv=c.moves; mv->i!=0; mv++) {
int i,j;
if (mv->when>0) continue;
i = mv->i; j = mv->j;
pt = (short *) vrna_alloc(sizeof(short)*(len+1));
memcpy(pt, c.pt,(len+1)*sizeof(short));
if (j<0) { /*it's a delete move */
pt[-i]=0;
pt[-j]=0;
} else { /* insert move */
if ((loopidx[i] == loopidx[j]) && /* i and j belong to same loop */
(pt[i] == 0) && (pt[j]==0) /* ... and are unpaired */
) {
pt[i]=j;
pt[j]=i;
} else {
free(pt);
continue; /* llegal move, try next; */
}
}
#ifdef LOOP_EN
en = c.curr_en + vrna_eval_move_pt(vc, c.pt, i, j);
#else
en = vrna_eval_structure_pt(vc, pt);
#endif
if (en<maxE) {
next[num_next].Sen = (en>oldE)?en:oldE;
next[num_next].curr_en = en;
next[num_next].pt = pt;
mv->when=dist;
mv->E = en;
next[num_next++].moves = copy_moves(c.moves);
mv->when=0;
}
else free(pt);
}
free(loopidx);
return num_next;
}
PRIVATE int find_path_once(vrna_fold_compound_t *vc, const char *struc1, const char *struc2, int maxE, int maxl) {
short *pt1, *pt2;
move_t *mlist;
int i, len, d, dist=0, result;
intermediate_t *current, *next;
pt1 = vrna_ptable(struc1);
pt2 = vrna_ptable(struc2);
len = (int) strlen(struc1);
mlist = (move_t *) vrna_alloc(sizeof(move_t)*len); /* bp_dist < n */
for (i=1; i<=len; i++) {
if (pt1[i] != pt2[i]) {
if (i<pt1[i]) { /* need to delete this pair */
mlist[dist].i = -i;
mlist[dist].j = -pt1[i];
mlist[dist++].when = 0;
}
if (i<pt2[i]) { /* need to insert this pair */
mlist[dist].i = i;
mlist[dist].j = pt2[i];
mlist[dist++].when = 0;
}
}
}
free(pt2);
BP_dist = dist;
current = (intermediate_t *) vrna_alloc(sizeof(intermediate_t)*(maxl+1));
current[0].pt = pt1;
current[0].Sen = current[0].curr_en = vrna_eval_structure_pt(vc, pt1);
current[0].moves = mlist;
next = (intermediate_t *) vrna_alloc(sizeof(intermediate_t)*(dist*maxl+1));
for (d=1; d<=dist; d++) { /* go through the distance classes */
int c, u, num_next=0;
intermediate_t *cc;
for (c=0; current[c].pt != NULL; c++) {
num_next += try_moves(vc, current[c], maxE, next+num_next, d);
}
if (num_next==0) {
for (cc=current; cc->pt != NULL; cc++) free_intermediate(cc);
current[0].Sen=INT_MAX;
break;
}
/* remove duplicates via sort|uniq
if this becomes a bottleneck we can use a hash instead */
qsort(next, num_next, sizeof(intermediate_t),compare_ptable);
for (u=0,c=1; c<num_next; c++) {
if (memcmp(next[u].pt,next[c].pt,sizeof(short)*len)!=0) {
next[++u] = next[c];
} else {
free_intermediate(next+c);
}
}
num_next = u+1;
qsort(next, num_next, sizeof(intermediate_t),compare_energy);
/* free the old stuff */
for (cc=current; cc->pt != NULL; cc++) free_intermediate(cc);
for (u=0; u<maxl && u<num_next; u++) {
current[u] = next[u];
}
for (; u<num_next; u++)
free_intermediate(next+u);
num_next=0;
}
free(next);
path = current[0].moves;
result = current[0].Sen;
free(current[0].pt); free(current);
return(result);
}
PRIVATE void free_intermediate(intermediate_t *i) {
free(i->pt);
free(i->moves);
i->pt = NULL;
i->moves = NULL;
i->Sen = INT_MAX;
}
PRIVATE int compare_ptable(const void *A, const void *B) {
intermediate_t *a, *b;
int c;
a = (intermediate_t *) A;
b = (intermediate_t *) B;
c = memcmp(a->pt, b->pt, a->pt[0]*sizeof(short));
if (c!=0) return c;
if ((a->Sen - b->Sen) != 0) return (a->Sen - b->Sen);
return (a->curr_en - b->curr_en);
}
PRIVATE int compare_energy(const void *A, const void *B) {
intermediate_t *a, *b;
a = (intermediate_t *) A;
b = (intermediate_t *) B;
if ((a->Sen - b->Sen) != 0) return (a->Sen - b->Sen);
return (a->curr_en - b->curr_en);
}
PRIVATE int compare_moves_when(const void *A, const void *B) {
move_t *a, *b;
a = (move_t *) A;
b = (move_t *) B;
return(a->when - b->when);
}
PRIVATE move_t* copy_moves(move_t *mvs) {
move_t *new;
new = (move_t *) vrna_alloc(sizeof(move_t)*(BP_dist+1));
memcpy(new,mvs,sizeof(move_t)*(BP_dist+1));
return new;
}
#ifdef TEST_FINDPATH
PUBLIC void print_path(const char *seq, const char *struc) {
int d;
char *s;
s = strdup(struc);
if (cut_point == -1)
printf("%s\n%s\n", seq, s);
/* printf("%s\n%s %6.2f\n", seq, s, vrna_eval_structure_simple(seq,s)); */
else {
char *pstruct, *pseq;
pstruct = vrna_cut_point_insert(s, cut_point);
pseq = vrna_cut_point_insert(seq, cut_point);
printf("%s\n%s\n", pseq, pstruct);
/* printf("%s\n%s %6.2f\n", pseq, pstruct, vrna_eval_structure_simple(seq,s)); */
free(pstruct);
free(pseq);
}
qsort(path, BP_dist, sizeof(move_t), compare_moves_when);
for (d=0; d<BP_dist; d++) {
int i,j;
i = path[d].i; j=path[d].j;
if (i<0) { /* delete */
s[(-i)-1] = s[(-j)-1] = '.';
} else {
s[i-1] = '('; s[j-1] = ')';
}
/* printf("%s %6.2f - %6.2f\n", s, vrna_eval_structure_simple(seq,s), path[d].E/100.0); */
}
free(s);
}
int main(int argc, char *argv[]) {
char *line, *seq, *s1, *s2;
int E, maxkeep=1000;
int verbose=0, i;
vrna_path_t *route, *r;
for (i=1; i<argc; i++) {
switch ( argv[i][1] ) {
case 'm': if (strcmp(argv[i],"-m")==0)
sscanf(argv[++i], "%d", &maxkeep);
break;
case 'v': verbose = !strcmp(argv[i],"-v");
break;
case 'd': if (strcmp(argv[i],"-d")==0){
sscanf(argv[++i], "%d", &dangles);
md.dangles = dangles;
}
break;
default: usage();
}
}
cut_point = -1;
line = vrna_read_line(stdin);
seq = vrna_cut_point_remove(line, &cut_point);
free(line);
line = vrna_read_line(stdin);
s1 = vrna_cut_point_remove(line, &cut_point);
free(line);
line = vrna_read_line(stdin);
s2 = vrna_cut_point_remove(line, &cut_point);
free(line);
E = find_saddle(seq, s1, s2, maxkeep);
printf("saddle_energy = %6.2f\n", E/100.);
if (verbose) {
if (path_fwd)
print_path(seq,s1);
else
print_path(seq,s2);
free(path);
path = NULL;
route = get_path(seq, s1, s2, maxkeep);
for (r=route; r->s; r++) {
if (cut_point == -1) {
printf("%s %6.2f\n", r->s, r->en);
/* printf("%s %6.2f - %6.2f\n", r->s, vrna_eval_structure_simple(seq,r->s), r->en); */
} else {
char *pstruct;
pstruct = vrna_cut_point_insert(r->s, cut_point);
printf("%s %6.2f\n", pstruct, r->en);
/* printf("%s %6.2f - %6.2f\n", pstruct, vrna_eval_structure_simple(seq,r->s), r->en); */
free(pstruct);
}
free(r->s);
}
free(route);
}
free(seq); free(s1); free(s2);
return(EXIT_SUCCESS);
}
static void usage(void){
vrna_message_error("usage: findpath.c [-m depth] [-d[0|1|2]] [-v]");
}
#endif