--- /dev/null
+/*
+ Last changed Time-stamp: <2010-06-24 15:12:09 ivo>
+ c Christoph Flamm and Ivo L Hofacker
+ {xtof,ivo}@tbi.univie.ac.at
+ Kinfold: $Name: $
+ $Id: baum.c,v 1.9 2008/05/21 10:15:45 ivo Exp $
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include "fold_vars.h"
+#include "fold.h"
+#include "utils.h"
+#include "pair_mat.h"
+#include "nachbar.h"
+#include "globals.h"
+
+#define MYTURN 1
+#define SAME_STRAND(I,J) (((I)>=cut_point)||((J)<cut_point))
+#define ORDER(x,y) if ((x)->nummer>(y)->nummer) {tempb=x; x=y; y=tempb;}
+
+/* item of structure ringlist */
+typedef struct _baum {
+ int nummer; /* number of base in sequence */
+ char typ; /* 'r' virtualroot, 'p' or 'q' paired, 'u' unpaired */
+ unsigned short base; /* 0<->unknown, 1<->A, 2<->C, 3<->G, 4<->U */
+ int loop_energy;
+ struct _baum *up;
+ struct _baum *next;
+ struct _baum *prev;
+ struct _baum *down;
+} baum;
+
+static char UNUSED rcsid[]="$Id: baum.c,v 1.9 2008/05/21 10:15:45 ivo Exp $";
+static short *pairList = NULL;
+static short *typeList = NULL;
+static short *aliasList = NULL;
+static baum *rl = NULL; /* ringlist */
+static baum *wurzl = NULL; /* virtualroot of ringlist-tree */
+static char **ptype = NULL;
+
+static int comp_struc(const void *A, const void *B);
+/* PUBLIC FUNCTIONES */
+void ini_or_reset_rl (void);
+void move_it (void);
+void update_tree (int i, int j);
+void clean_up_rl (void);
+
+/* PRIVATE FUNCTIONES */
+static void ini_ringlist(void);
+static void reset_ringlist(void);
+static void struc2tree (char *struc);
+static void close_bp_en (baum *i, baum *j);
+static void close_bp (baum *i, baum *j);
+static void open_bp (baum *i);
+static void open_bp_en (baum *i);
+static void inb (baum *root);
+static void inb_nolp (baum *root);
+static void dnb (baum *rli);
+static void dnb_nolp (baum *rli);
+static void fnb (baum *rli);
+static void make_ptypes(const short *S);
+/* debugging tool(s) */
+#if 0
+static void rl_status(void);
+#endif
+
+/* convert structure in bracked-dot-notation to a ringlist-tree */
+static void struc2tree(char *struc) {
+ char* struc_copy;
+ int ipos, jpos, balance = 0;
+ baum *rli, *rlj;
+
+ struc_copy = (char *)calloc(GSV.len+1, sizeof(char));
+ assert(struc_copy);
+ strcpy(struc_copy,struc);
+
+ for (ipos = 0; ipos < GSV.len; ipos++) {
+ if (struc_copy[ipos] == ')') {
+ jpos = ipos;
+ struc_copy[ipos] = '.';
+ balance++;
+ while (struc_copy[--ipos] != '(');
+ struc_copy[ipos] = '.';
+ balance--;
+ rli = &rl[ipos];
+ rlj = &rl[jpos];
+ close_bp(rli, rlj);
+ }
+ }
+
+ if (balance) {
+ fprintf(stderr,
+ "struc2tree(): start structure is not balanced !\n%s\n%s\n",
+ GAV.farbe, struc);
+ exit(1);
+ }
+
+ GSV.currE = GSV.startE =
+ (float )energy_of_structure_pt(GAV.farbe,
+ pairList, typeList, aliasList, 0) / 100.0;
+ {
+ int i;
+ for(i = 0; i < GSV.len; i++) {
+ if (pairList[i+1]>i+1)
+ rl[i].loop_energy = loop_energy(pairList, typeList, aliasList,i+1);
+ }
+ wurzl->loop_energy = loop_energy(pairList, typeList, aliasList,0);
+ }
+
+ free(struc_copy);
+}
+
+/**/
+static void ini_ringlist(void) {
+ int i;
+
+ /* needed by function energy_of_struct_pt() from Vienna-RNA-1.4 */
+ pairList = (short *)calloc(GSV.len + 2, sizeof(short));
+ assert(pairList != NULL);
+ typeList = (short *)calloc(GSV.len + 2, sizeof(short));
+ assert(typeList != NULL);
+ aliasList = (short *)calloc(GSV.len + 2, sizeof(short));
+ assert(aliasList != NULL);
+ pairList[0] = typeList[0] = aliasList[0] = GSV.len;
+ ptype = (char **)calloc(GSV.len + 2, sizeof(char *));
+ assert(ptype != NULL);
+ for (i=0; i<=GSV.len; i++) {
+ ptype[i] = (char*)calloc(GSV.len + 2, sizeof(char));
+ assert(ptype[i] != NULL);
+ }
+
+ /* allocate virtual root */
+ wurzl = (baum *)calloc(1, sizeof(baum));
+ assert(wurzl != NULL);
+ /* allocate ringList */
+ rl = (baum *)calloc(GSV.len+1, sizeof(baum));
+ assert(rl != NULL);
+ /* allocate PostOrderList */
+
+ /* initialize virtualroot */
+ wurzl->typ = 'r';
+ wurzl->nummer = -1;
+ /* connect virtualroot to ringlist-tree in down direction */
+ wurzl->down = &rl[GSV.len];
+ /* initialize post-order list */
+
+ make_pair_matrix();
+
+ /* initialize rest of ringlist-tree */
+ for(i = 0; i < GSV.len; i++) {
+ int c;
+ GAV.currform[i] = '.';
+ GAV.prevform[i] = 'x';
+ pairList[i+1] = 0;
+ rl[i].typ = 'u';
+ /* decode base to numeric value */
+ c = encode_char(GAV.farbe[i]);
+ rl[i].base = typeList[i+1] = c;
+ aliasList[i+1] = alias[typeList[i+1]];
+ /* astablish links for node of the ringlist-tree */
+ rl[i].nummer = i;
+ rl[i].next = &rl[i+1];
+ rl[i].prev = ((i == 0) ? &rl[GSV.len] : &rl[i-1]);
+ rl[i].up = rl[i].down = NULL;
+ }
+ GAV.currform[GSV.len] = GAV.prevform[GSV.len] = '\0';
+ make_ptypes(aliasList);
+
+ rl[i].nummer = i;
+ rl[i].base = 0;
+ /* make ringlist circular in next, prev direction */
+ rl[i].next = &rl[0];
+ rl[i].prev = &rl[i-1];
+ /* make virtual basepair for virtualroot */
+ rl[i].up = wurzl;
+ rl[i].typ = 'x';
+
+}
+
+/**/
+void ini_or_reset_rl(void) {
+
+ /* if there is no ringList-tree make a new one */
+ if (wurzl == NULL) {
+ ini_ringlist();
+
+ /* start structure */
+ struc2tree(GAV.startform);
+ GSV.currE = GSV.startE = energy_of_structure(GAV.farbe, GAV.startform, 0);
+
+
+ /* stop structure(s) */
+ if ( GTV.stop ) {
+ int i;
+
+ qsort(GAV.stopform, GSV.maxS, sizeof(char *), comp_struc);
+ for (i = 0; i< GSV.maxS; i++)
+ GAV.sE[i] = energy_of_structure(GAV.farbe_full, GAV.stopform[i], 0);
+ }
+ else {
+ if(GTV.noLP)
+ noLonelyPairs=1;
+ initialize_cofold(GSV.len);
+ /* fold sequence to get Minimum free energy structure (Mfe) */
+ GAV.sE[0] = cofold(GAV.farbe_full, GAV.stopform[0]);
+ free_arrays();
+ /* revaluate energy of Mfe (maye differ if --logML=logarthmic */
+ GAV.sE[0] = energy_of_structure(GAV.farbe_full, GAV.stopform[0], 0);
+ }
+ GSV.stopE = GAV.sE[0];
+ ini_nbList(strlen(GAV.farbe_full)*strlen(GAV.farbe_full));
+ }
+ else {
+ /* reset ringlist-tree to start conditions */
+ reset_ringlist();
+ if(GTV.start) struc2tree(GAV.startform);
+ else {
+ GSV.currE = GSV.startE;
+ }
+ }
+}
+
+/**/
+static void reset_ringlist(void) {
+ int i;
+
+ for(i = 0; i < GSV.len; i++) {
+ GAV.currform[i] = '.';
+ GAV.prevform[i] = 'x';
+ pairList[i+1] = 0;
+ rl[i].typ = 'u';
+ rl[i].next = &rl[i + 1];
+ rl[i].prev = ((i == 0) ? &rl[GSV.len] : &rl[i - 1]);
+ rl[i].up = rl[i].down = NULL;
+ }
+ rl[i].next = &rl[0];
+ rl[i].prev = &rl[i-1];
+ rl[i].up = wurzl;
+}
+
+/* update ringlist-tree */
+void update_tree(int i, int j) {
+
+ baum *rli, *rlj, *tempb;
+
+ if ( abs(i) < GSV.len) { /* >> single basepair move */
+ if ((i > 0) && (j > 0)) { /* insert */
+ rli = &rl[i-1];
+ rlj = &rl[j-1];
+ close_bp_en(rli, rlj);
+ }
+ else if ((i < 0)&&(j < 0)) { /* delete */
+ i = -i;
+ rli = &rl[i-1];
+ open_bp_en(rli);
+ }
+ else { /* shift */
+ if (i > 0) { /* i remains the same, j shifts */
+ j=-j;
+ rli=&rl[i-1];
+ rlj=&rl[j-1];
+ open_bp_en(rli);
+ ORDER(rli, rlj);
+ close_bp_en(rli, rlj);
+ }
+ else { /* j remains the same, i shifts */
+ baum *old_rli;
+ i = -i;
+ rli = &rl[i-1];
+ rlj = &rl[j-1];
+ old_rli = rlj->up;
+ open_bp_en(old_rli);
+ ORDER(rli, rlj);
+ close_bp_en(rli, rlj);
+ }
+ }
+ } /* << single basepair move */
+ else { /* >> double basepair move */
+ if ((i > 0) && (j > 0)) { /* insert */
+ rli = &rl[i-GSV.len-2];
+ rlj = &rl[j-GSV.len-2];
+ close_bp_en(rli->next, rlj->prev);
+ close_bp_en(rli, rlj);
+ }
+ else if ((i < 0)&&(j < 0)) { /* delete */
+ i = -i;
+ rli = &rl[i-GSV.len-2];
+ open_bp_en(rli);
+ open_bp_en(rli->next);
+ }
+ } /* << double basepair move */
+
+}
+
+/* open a particular base pair */
+void open_bp(baum *i) {
+
+ baum *in; /* points to i->next */
+
+ /* change string representation */
+ GAV.currform[i->nummer] = '.';
+ GAV.currform[i->down->nummer] = '.';
+
+ /* change pairtable representation */
+ pairList[1 + i->nummer] = 0;
+ pairList[1 + i->down->nummer] = 0;
+
+ /* change tree representation */
+ in = i->next;
+ i->typ = 'u';
+ i->down->typ = 'u';
+ i->next = i->down->next;
+ i->next->prev = i;
+ in->prev = i->down;
+ i->down->next = in;
+ i->down = in->prev->up = NULL;
+}
+
+/* close a particular base pair */
+void close_bp (baum *i, baum *j) {
+
+ baum *jn; /* points to j->next */
+
+ /* change string representation */
+ GAV.currform[i->nummer] = '(';
+ GAV.currform[j->nummer] = ')';
+
+ /* change pairtable representation */
+ pairList[1 + i->nummer] = 1+ j->nummer;
+ pairList[1 + j->nummer] = 1 + i->nummer;
+
+ /* change tree representation */
+ jn = j->next;
+ i->typ = 'p';
+ j->typ = 'q';
+ i->down = j;
+ j->up = i;
+ i->next->prev = j;
+ j->next->prev = i;
+ j->next = i->next;
+ i->next = jn;
+}
+
+# if 0
+/* for a given tree, generate postorder-list */
+static void make_poList (baum *root) {
+
+ baum *stop, *rli;
+
+ if (!root) root = wurzl;
+ stop = root->down;
+
+ /* foreach base in ringlist ... */
+ for (rli = stop->next; rli != stop; rli = rli->next) {
+ /* ... explore subtee if bp found */
+ if (rli->typ == 'p') {
+ /* fprintf(stderr, "%d >%d<\n", poListop, rli->nummer); */
+ poList[poListop++] = rli;
+ if ( poListop > GSV.len+1 ) {
+ fprintf(stderr, "Something went wrong in make_poList()\n");
+ exit(1);
+ }
+ make_poList(rli);
+ }
+ }
+ return;
+}
+#endif
+
+/* for a given ringlist, generate all structures
+ with one additional basepair */
+static void inb(baum *root) {
+
+ int EoT;
+ int E_old, E_new_in, E_new_out;
+ baum *stop,*rli,*rlj;
+
+ E_old = root->loop_energy;
+ stop=root->down;
+ /* loop ringlist over all possible i positions */
+ for(rli=stop->next;rli!=stop;rli=rli->next){
+ /* potential i-position is already paired */
+ if(rli->typ=='p') continue;
+ /* loop ringlist over all possible j positions */
+ for(rlj=rli->next;rlj!=stop;rlj=rlj->next){
+ /* base pair must enclose at least 3 bases */
+ if(rlj->nummer - rli->nummer < MYTURN) continue;
+ /* potential j-position is already paired */
+ if(rlj->typ=='p') continue;
+ /* if i-j can form a base pair ... */
+ if(ptype[rli->nummer][rlj->nummer]){
+ /* close the base bair and ... */
+ close_bp(rli,rlj);
+ E_new_in = loop_energy(pairList, typeList, aliasList,rli->nummer+1);
+ E_new_out = loop_energy(pairList, typeList, aliasList,root->nummer+1);
+ /* ... evaluate energy of the structure */
+ EoT = (int) (GSV.currE*100 + ((GSV.currE<0)?-0.4:0.4)) + E_new_in + E_new_out - E_old ;
+ /* assert(EoT == energy_of_struct_pt(GAV.farbe, pairList, typeList, aliasList)); */
+ /* open the base pair again... */
+ open_bp(rli);
+ /* ... and put the move and the enegy
+ of the structure into the neighbour list */
+ update_nbList(1 + rli->nummer, 1 + rlj->nummer, EoT);
+ }
+ }
+ }
+}
+
+/* for a given ringlist, generate all structures (canonical)
+ with one additional base pair (BUT WITHOUT ISOLATED BASE PAIRS) */
+static void inb_nolp(baum *root) {
+
+ int EoT = 0;
+ baum *stop, *rli, *rlj;
+
+ stop = root->down;
+ /* loop ringlist over all possible i positions */
+ for (rli=stop->next;rli!=stop;rli=rli->next) {
+ /* potential i-position is already paired */
+ if (rli->typ=='p') continue;
+ /* loop ringlist over all possible j positions */
+ for (rlj=rli->next;rlj!=stop;rlj=rlj->next) {
+ /* base pair must enclose at least 3 bases */
+ if (rlj->nummer - rli->nummer < MYTURN) continue;
+ /* potential j-position is already paired */
+ if (rlj->typ=='p') continue;
+ /* if i-j can form a base pair ... */
+ if (ptype[rli->nummer][rlj->nummer]) {
+ /* ... and extends a helix ... */
+ if (((rli->prev==stop && rlj->next==stop) && stop->typ != 'x') ||
+ (rli->next == rlj->prev)) {
+ /* ... close the base bair and ... */
+ close_bp(rli,rlj);
+ /* ... evaluate energy of the structure */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* open the base pair again... */
+ open_bp(rli);
+ /* ... and put the move and the enegy
+ of the structure into the neighbour list */
+ update_nbList(1 + rli->nummer, 1 + rlj->nummer, EoT);
+ }
+ /* if double insertion is possible ... */
+ else if ((rlj->nummer - rli->nummer >= MYTURN+2)&&
+ (rli->next->typ != 'p' && rlj->prev->typ != 'p') &&
+ (rli->next->next != rlj->prev->prev) &&
+ (ptype[rli->next->nummer][rlj->prev->nummer])) {
+ /* close the two base bair and ... */
+ close_bp(rli->next, rlj->prev);
+ close_bp(rli, rlj);
+ /* ... evaluate energy of the structure */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* open the two base pair again ... */
+ open_bp(rli);
+ open_bp(rli->next);
+ /* ... and put the move and the enegy
+ of the structure into the neighbour list */
+ update_nbList(1+rli->nummer+GSV.len+1, 1+rlj->nummer+GSV.len+1, EoT);
+ }
+ }
+ }
+ }
+}
+
+/* for a given ringlist, generate all structures
+ with one less base pair */
+static void dnb(baum *rli){
+
+ int EoT, E_old_in, E_old_out, E_new;
+
+ baum *rlj, *r;
+
+ rlj=rli->down;
+ open_bp(rli);
+ /* ... evaluate energy of the structure */
+
+ for (r=rli->next; r->up==NULL; r=r->next);
+ E_old_in = rli->loop_energy;
+ E_old_out = r->up->loop_energy;
+ E_new = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
+ EoT = (int) (GSV.currE*100 + ((GSV.currE<0)?-0.4:0.4)) -
+ E_old_in - E_old_out + E_new;
+
+ /* assert(EoT== energy_of_struct_pt(GAV.farbe, pairList, typeList, aliasList));*/
+ close_bp(rli,rlj);
+ update_nbList(-(1 + rli->nummer), -(1 + rlj->nummer), EoT);
+}
+
+/* for a given ringlist, generate all structures (canonical)
+ with one less base pair (BUT WITHOUT ISOLATED BASE PAIRS) */
+static void dnb_nolp(baum *rli) {
+
+ int EoT = 0;
+ baum *rlj;
+ baum *rlin = NULL; /* pointers to following pair in helix, if any */
+ baum *rljn = NULL;
+ baum *rlip = NULL; /* pointers to preceding pair in helix, if any */
+ baum *rljp = NULL;
+
+ rlj = rli->down;
+
+ /* immediate interior base pair ? */
+ if (rlj->next == rlj->prev) {
+ rlin = rlj->next;
+ rljn = rlin->down;
+ }
+
+ /* immediate exterior base pair and not virtualroot ? */
+ if (rli->prev == rli->next && rli->next->typ != 'x') {
+ rlip = rli->next->up;
+ rljp = rli->next;
+ }
+
+ /* double delete ? */
+ if (rlip==NULL && rlin && rljn->next != rljn->prev ) {
+ /* open the two base pairs ... */
+ open_bp(rli);
+ open_bp(rlin);
+ /* ... evaluate energy of the structure ... */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* ... and put the move and the enegy
+ of the structure into the neighbour list ... */
+ update_nbList(-(1+rli->nummer+GSV.len+1),-(1+rlj->nummer+GSV.len+1), EoT);
+ /* ... and close the two base pairs again */
+ close_bp(rlin, rljn);
+ close_bp(rli, rlj);
+ } else { /* single delete */
+ /* the following will work only if boolean expr are shortcicuited */
+ if (rlip==NULL || (rlip->prev == rlip->next && rlip->prev->typ != 'x'))
+ if (rlin ==NULL || (rljn->next == rljn->prev)) {
+ /* open the base pair ... */
+ open_bp(rli);
+ /* ... evaluate energy of the structure ... */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* ... and put the move and the enegy
+ of the structure into the neighbour list ... */
+ update_nbList(-(1 + rli->nummer),-(1 + rlj->nummer), EoT);
+ /* and close the base pair again */
+ close_bp(rli, rlj);
+ }
+ }
+}
+
+/* for a given ringlist, generate all structures
+ with one shifted base pair */
+static void fnb(baum *rli) {
+
+ int EoT = 0, x;
+ baum *rlj, *stop, *help_rli, *help_rlj;
+
+ stop = rli->down;
+
+ /* examin interior loop of bp(ij); (.......)
+ i of j move -> <- */
+ for (rlj = stop->next; rlj != stop; rlj = rlj->next) {
+ /* prevent shifting to paired position */
+ if ((rlj->typ=='p')||(rlj->typ=='q')) continue;
+ /* j-position of base pair shifts to k position (ij)->(ik) i<k<j */
+ if ( (rlj->nummer-rli->nummer >= MYTURN)
+ && (ptype[rli->nummer][rlj->nummer]) ) {
+ /* open original basepair */
+ open_bp(rli);
+ /* close shifted version of original basepair */
+ close_bp(rli, rlj);
+ /* evaluate energy of the structure */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* put the move and the enegy of the structure into the neighbour list */
+ update_nbList(1+rli->nummer, -(1+rlj->nummer), EoT);
+ /* open shifted basepair */
+ open_bp(rli);
+ /* restore original basepair */
+ close_bp(rli, stop);
+ }
+ /* i-position of base pair shifts to position k (ij)->(kj) i<k<j */
+ if ( (stop->nummer-rlj->nummer >= MYTURN)
+ && (ptype[stop->nummer][rlj->nummer]) ) {
+ /* open original basepair */
+ open_bp(rli);
+ /* close shifted version of original basepair */
+ close_bp(rlj, stop);
+ /* evaluate energy of the structure */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* put the move and the enegy of the structure into the neighbour list */
+ update_nbList(-(1 + rlj->nummer), 1 + stop->nummer, EoT);
+ /* open shifted basepair */
+ open_bp(rlj);
+ /* restore original basepair */
+ close_bp(rli, stop);
+ }
+ }
+ /* examin exterior loop of bp(ij); (.......)
+ i or j moves <- -> */
+ for (rlj=rli->next;rlj!=rli;rlj=rlj->next) {
+ if ((rlj->typ=='p') || (rlj->typ=='q' ) || (rlj->typ=='x')) continue;
+ x=rlj->nummer-rli->nummer;
+ if (x<0) x=-x;
+ /* j-position of base pair shifts to position k */
+ if ((x >= MYTURN) && (ptype[rli->nummer][rlj->nummer])) {
+ if (rli->nummer<rlj->nummer) {
+ help_rli=rli;
+ help_rlj=rlj;
+ }
+ else {
+ help_rli=rlj;
+ help_rlj=rli;
+ }
+ /* open original basepair */
+ open_bp(rli);
+ /* close shifted version of original basepair */
+ close_bp(help_rli,help_rlj);
+ /* evaluate energy of the structure */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* put the move and the enegy of the structure into the neighbour list */
+ update_nbList(1 + rli->nummer, -(1 + rlj->nummer), EoT);
+ /* open shifted base pair */
+ open_bp(help_rli);
+ /* restore original basepair */
+ close_bp(rli,stop);
+ }
+ x = rlj->nummer-stop->nummer;
+ if (x < 0) x = -x;
+ /* i-position of base pair shifts to position k */
+ if ((x >= MYTURN) && (ptype[stop->nummer][rlj->nummer])) {
+ if (stop->nummer < rlj->nummer) {
+ help_rli = stop;
+ help_rlj = rlj;
+ }
+ else {
+ help_rli = rlj;
+ help_rlj = stop;
+ }
+ /* open original basepair */
+ open_bp(rli);
+ /* close shifted version of original basepair */
+ close_bp(help_rli, help_rlj);
+ /* evaluate energy of the structure */
+ EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
+ /* put the move and the enegy of the structure into the neighbour list */
+ update_nbList(-(1 + rlj->nummer), 1 + stop->nummer, EoT);
+ /* open shifted basepair */
+ open_bp(help_rli);
+ /* restore original basepair */
+ close_bp(rli,stop);
+ }
+ }
+}
+
+/* for a given tree (structure),
+ generate all neighbours according to moveset */
+void move_it (void) {
+ int i;
+
+ GSV.currE =
+ energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0)/100.;
+
+ if ( GTV.noLP ) { /* canonical neighbours only */
+ inb_nolp(wurzl);
+ for (i = 0; i < GSV.len; i++) {
+
+ if (pairList[i+1]>i+1) {
+ inb_nolp(rl+i); /* insert pair neighbours */
+ dnb_nolp(rl+i); /* delete pair neighbour */
+ }
+ }
+ }
+ else { /* all neighbours */
+ inb(wurzl);
+ for (i = 0; i < GSV.len; i++) {
+
+ if (pairList[i+1]>i+1) {
+ inb(rl+i); /* insert pair neighbours */
+ dnb(rl+i); /* delete pair neighbour */
+ if ( GTV.noShift == 0 ) fnb(rl+i);
+ }
+ }
+ }
+}
+
+
+/**/
+void clean_up_rl(void) {
+ int i;
+ free(pairList); pairList=NULL;
+ free(typeList); typeList = NULL;
+ free(aliasList); aliasList = NULL;
+ free(rl); rl=NULL;
+ free(wurzl); wurzl=NULL;
+ for (i=0; i<=GSV.len; i++)
+ free(ptype[i]);
+ free(ptype);
+ ptype=NULL;
+}
+
+/**/
+static int comp_struc(const void *A, const void *B) {
+ int aE, bE;
+ aE = (int)(100 * energy_of_structure(GAV.farbe_full, ((char **)A)[0], 0));
+ bE = (int)(100 * energy_of_structure(GAV.farbe_full, ((char **)B)[0], 0));
+ return (aE-bE);
+}
+
+#if 0
+/**/
+static void rl_status(void) {
+
+ int i;
+
+ printf("\n%s\n%s\n", GAV.farbe, GAV.currform);
+ for (i=0; i <= GSV.len; i++) {
+ printf("%2d %c %c %2d %2d %2d %2d\n",
+ rl[i].nummer,
+ i == GSV.len ? 'X': GAV.farbe[i],
+ rl[i].typ,
+ rl[i].up==NULL?0:(rl[i].up)->nummer,
+ rl[i].down==NULL?0:(rl[i].down)->nummer,
+ (rl[i].prev)->nummer,
+ (rl[i].next)->nummer);
+ }
+ printf("---\n");
+}
+#endif
+
+#define TURN 3
+static void make_ptypes(const short *S) {
+ int n,i,j,k,l;
+ n=S[0];
+ for (k=1; k<n; k++)
+ for (l=1; l<=2; l++) {
+ int type,ntype=0,otype=0;
+ i=k; j = i+l+TURN;
+ if (!SAME_STRAND(i,j)) j=cut_point;
+ if (j>n) continue;
+ type = pair[S[i]][S[j]];
+ while ((i>=1)&&(j<=n)) {
+ if ((i>1)&&(j<n)) ntype = pair[S[i-1]][S[j+1]];
+ if (noLonelyPairs && (!otype) && (!ntype))
+ type = 0; /* i.j can only form isolated pairs */
+ ptype[i-1][j-1] = ptype[j-1][i-1] = (char) type;
+ otype = type;
+ type = ntype;
+ i--; j++;
+ }
+ }
+}
+
+static void close_bp_en (baum *i, baum *j) {
+ /* close bp and update energy */
+ baum *r;
+ close_bp(i,j);
+ i->loop_energy = loop_energy(pairList,typeList,aliasList,i->nummer+1);
+ for (r=i->next; r->up==NULL; r=r->next);
+ r->up->loop_energy = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
+};
+
+static void open_bp_en (baum *i) {
+ /* open bp and update energy */
+ baum *r;
+ i->loop_energy=0;
+ open_bp(i);
+ for (r=i->next; r->up==NULL; r=r->next);
+ r->up->loop_energy = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
+};
git://source.jalview.org / jabaws.git / blobdiff