--- /dev/null
+
+/* Last changed Time-stamp: <2007-12-05 14:05:51 ivo> */
+/*
+ minimum free energy
+ RNA secondary structure prediction
+
+ c Ivo Hofacker, Chrisoph Flamm
+ original implementation by
+ Walter Fontana
+
+ Vienna RNA package
+*/
+
+#include <config.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <ctype.h>
+#include <string.h>
+#include "utils.h"
+#include "energy_par.h"
+#include "fold_vars.h"
+#include "pair_mat.h"
+#include "params.h"
+#include "snofold.h"
+#include "loop_energies.h"
+/*@unused@*/
+static char rcsid[] UNUSED = "$Id: fold.c,v 1.38 2007/12/19 10:27:42 ivo Exp $";
+#ifdef __GNUC__
+#define INLINE inline
+#else
+#define INLINE
+#endif
+
+#define PAREN
+
+#define PUBLIC
+#define PRIVATE static
+
+#define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
+#define NEW_NINIO 1 /* new asymetry penalty */
+
+/*@unused@*/
+PRIVATE void letter_structure(char *structure, int length) UNUSED;
+PRIVATE void parenthesis_structure(char *structure, int length);
+PRIVATE void get_arrays(unsigned int size);
+/* PRIVATE void scale_parameters(void); */
+/* PRIVATE int stack_energy(int i, const char *string); */
+PRIVATE void make_ptypes(const short *S, const char *structure);
+PRIVATE void encode_seq(const char *sequence);
+PRIVATE void backtrack(const char *sequence, int s);
+PRIVATE int fill_arrays(const char *sequence, const int max_asymm, const int threshloop,
+ const int min_s2, const int max_s2, const int half_stem, const int max_half_stem);
+/*@unused@*/
+
+
+/* alifold */
+PRIVATE void alisnoinitialize_fold(const int length);
+PRIVATE void make_pscores(const short *const* S, const char *const* AS,int n_seq, const char *structure);
+PRIVATE int *pscore; /* precomputed array of pair types */
+PRIVATE short **Sali;
+PRIVATE int alifill_arrays(const char **string, const int max_asymm, const int threshloop,
+ const int min_s2, const int max_s2, const int half_stem,
+ const int max_half_stem);
+PRIVATE void aliget_arrays(unsigned int size);
+PRIVATE short * aliencode_seq(const char *sequence);
+PRIVATE int alibacktrack(const char **strings, int s);
+
+#define UNIT 100
+#define MINPSCORE -2 * UNIT
+/* end alifold */
+
+#define MAXSECTORS 500 /* dimension for a backtrack array */
+#define LOCALITY 0. /* locality parameter for base-pairs */
+
+#define MIN2(A, B) ((A) < (B) ? (A) : (B))
+#define MAX2(A, B) ((A) > (B) ? (A) : (B))
+#define SAME_STRAND(I,J) (((I)>=cut_point)||((J)<cut_point))
+
+PRIVATE paramT *P = NULL;
+
+PRIVATE int *indx = NULL; /* index for moving in the triangle matrices c[] and fMl[]*/
+
+PRIVATE int *c = NULL; /* energy array, given that i-j pair */
+PRIVATE int *cc = NULL; /* linear array for calculating canonical structures */
+PRIVATE int *cc1 = NULL; /* " " */
+PRIVATE int *Fmi = NULL; /* holds row i of fML (avoids jumps in memory) */
+PRIVATE int *DMLi = NULL; /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j]) */
+PRIVATE int *DMLi1 = NULL; /* MIN(fML[i+1,k]+fML[k+1,j]) */
+PRIVATE int *DMLi2 = NULL; /* MIN(fML[i+2,k]+fML[k+1,j]) */
+PRIVATE char *ptype = NULL; /* precomputed array of pair types */
+PRIVATE short *S = NULL, *S1 = NULL;
+PRIVATE int init_length=-1;
+PRIVATE int *mLoop = NULL; /*contains the minimum of c for a xy range*/
+PRIVATE folden **foldlist = NULL;
+PRIVATE folden **foldlist_XS = NULL;
+
+PRIVATE int *BP = NULL; /* contains the structure constrainsts: BP[i]
+ -1: | = base must be paired
+ -2: < = base must be paired with j<i
+ -3: > = base must be paired with j>i
+ -4: x = base must not pair
+ positive int: base is paired with int */
+
+
+static sect sector[MAXSECTORS]; /* stack of partial structures for backtracking */
+
+PRIVATE char alpha[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
+/* needed by cofold/eval */
+/* PRIVATE int cut_in_loop(int i); */
+/* PRIVATE int min_hairpin = TURN; */
+
+/* some definitions to take circfold into account... */
+/* PRIVATE int *fM2 = NULL;*/ /* fM2 = multiloop region with exactly two stems, extending to 3' end */
+PUBLIC int Fc, FcH, FcI, FcM; /* parts of the exterior loop energies */
+/*--------------------------------------------------------------------------*/
+
+void snoinitialize_fold(const int length)
+{
+ unsigned int n;
+ if (length<1) nrerror("snoinitialize_fold: argument must be greater 0");
+ if (init_length>0) snofree_arrays(length);
+ get_arrays((unsigned) length);
+ init_length=length;
+ for (n = 1; n <= (unsigned) length; n++)
+ indx[n] = (n*(n-1)) >> 1; /* n(n-1)/2 */
+
+ snoupdate_fold_params();
+}
+
+PRIVATE void alisnoinitialize_fold(const int length)
+{
+ unsigned int n;
+ if (length<1) nrerror("snoinitialize_fold: argument must be greater 0");
+ if (init_length>0) snofree_arrays(length);
+ aliget_arrays((unsigned) length);
+ make_pair_matrix();
+ init_length=length;
+ for (n = 1; n <= (unsigned) length; n++)
+ indx[n] = (n*(n-1)) >> 1; /* n(n-1)/2 */
+ snoupdate_fold_params();
+}
+
+
+/*--------------------------------------------------------------------------*/
+
+PRIVATE void get_arrays(unsigned int size)
+{
+ indx = (int *) space(sizeof(int)*(size+1));
+ c = (int *) space(sizeof(int)*((size*(size+1))/2+2));
+ mLoop = (int *) space(sizeof(int)*((size*(size+1))/2+2));
+
+ ptype = (char *) space(sizeof(char)*((size*(size+1))/2+2));
+ cc = (int *) space(sizeof(int)*(size+2));
+ cc1 = (int *) space(sizeof(int)*(size+2));
+ Fmi = (int *) space(sizeof(int)*(size+1));
+ DMLi = (int *) space(sizeof(int)*(size+1));
+ DMLi1 = (int *) space(sizeof(int)*(size+1));
+ DMLi2 = (int *) space(sizeof(int)*(size+1));
+ if (base_pair) free(base_pair);
+ base_pair = (struct bondT *) space(sizeof(struct bondT)*(1+size/2));
+ /* extra array(s) for circfold() */
+}
+
+PRIVATE void aliget_arrays(unsigned int size)
+{
+ indx = (int *) space(sizeof(int)*(size+1));
+ c = (int *) space(sizeof(int)*((size*(size+1))/2+2));
+ mLoop = (int *) space(sizeof(int)*((size*(size+1))/2+2));
+ pscore = (int *) space(sizeof(int)*((size*(size+1))/2+2));
+ ptype = (char *) space(sizeof(char)*((size*(size+1))/2+2));
+ cc = (int *) space(sizeof(int)*(size+2));
+ cc1 = (int *) space(sizeof(int)*(size+2));
+ Fmi = (int *) space(sizeof(int)*(size+1));
+ DMLi = (int *) space(sizeof(int)*(size+1));
+ DMLi1 = (int *) space(sizeof(int)*(size+1));
+ DMLi2 = (int *) space(sizeof(int)*(size+1));
+ if (base_pair) free(base_pair);
+ base_pair = (struct bondT *) space(sizeof(struct bondT)*(1+size/2));
+ /* extra array(s) for circfold() */
+}
+
+
+
+
+/*--------------------------------------------------------------------------*/
+
+
+void snofree_arrays(const int length)
+{
+ free(indx); free(c);free(cc); free(cc1);
+ free(ptype);free(mLoop);
+ int i;
+ for(i=length;i>-1;i--){
+ while(foldlist[i]!=NULL){
+ folden *n = foldlist[i];
+ foldlist[i] = foldlist[i]->next;
+ free(n);
+ }
+ free(foldlist[i]);
+ }
+ free(foldlist);
+ for(i=length;i>-1;i--){
+ while(foldlist_XS[i]!=NULL){
+ folden *n = foldlist_XS[i];
+ foldlist_XS[i] = foldlist_XS[i]->next;
+ free(n);
+ }
+ free(foldlist_XS[i]);
+ }
+ free(foldlist_XS);
+ free(base_pair); base_pair=NULL; free(Fmi);
+ free(DMLi); free(DMLi1);free(DMLi2);
+ free(BP);
+ init_length=0;
+}
+
+void alisnofree_arrays(const int length)
+{
+ free(indx); free(c);free(cc); free(cc1);
+ free(ptype);free(mLoop);free(pscore);
+ int i;
+ for(i=length-1;i>-1;i--){
+ while(foldlist[i]!=NULL){
+ folden *n = foldlist[i];
+ foldlist[i] = foldlist[i]->next;
+ free(n);
+ }
+ free(foldlist[i]);
+ }
+ free(foldlist);
+ free(base_pair); base_pair=NULL; free(Fmi);
+ free(DMLi); free(DMLi1);free(DMLi2);
+ free(BP);
+ init_length=0;
+}
+
+/*--------------------------------------------------------------------------*/
+
+void snoexport_fold_arrays(int **indx_p, int **mLoop_p, int **cLoop, folden ***fold_p, folden ***fold_p_XS) {
+ /* make the DP arrays available to routines such as subopt() */
+ *indx_p = indx; *mLoop_p = mLoop;
+ *cLoop = c; *fold_p = foldlist;*fold_p_XS=foldlist_XS;
+}
+
+/* void alisnoexport_fold_arrays(int **indx_p, int **mLoop_p, int **cLoop, folden ***fold_p, int **pscores) { */
+/* /\* make the DP arrays available to routines such as subopt() *\/ */
+/* *indx_p = indx; *mLoop_p = mLoop; */
+/* *cLoop = c; *fold_p = foldlist; */
+/* *pscores=pscore; */
+/* } */
+
+/*--------------------------------------------------------------------------*/
+
+
+
+
+
+
+
+
+
+
+int snofold(const char *string, char *structure, const int max_assym, const int threshloop,
+ const int min_s2, const int max_s2, const int half_stem, const int max_half_stem) {
+ int length, energy, bonus, bonus_cnt, s;
+
+ /* Variable initialization */
+ bonus = 0;
+ bonus_cnt = 0;
+ s = 0;
+ length = (int) strlen(string);
+
+ S = encode_sequence(string, 0);
+ S1 = encode_sequence(string, 1);
+
+
+ /* structure = (char *) space((unsigned) length+1); */
+
+ if (length>init_length) snoinitialize_fold(length);
+ else if (fabs(P->temperature - temperature)>1e-6) snoupdate_fold_params();
+
+
+
+ /* encode_seq(string); */
+ BP = (int *)space(sizeof(int)*(length+2));
+ make_ptypes(S, structure);
+ energy=fill_arrays(string, max_assym, threshloop, min_s2, max_s2, half_stem, max_half_stem);
+ backtrack(string, s);
+
+#if 0
+ /*no structure output, no backtrack*/
+ backtrack(string, 0);
+ parenthesis_structure(structure, length);
+#endif
+ free(structure);
+ free(S); free(S1); /* free(BP); */
+ return energy;
+}
+
+PRIVATE void make_pscores(const short *const* S, const char *const* AS,
+ int n_seq, const char *structure) {
+ /* calculate co-variance bonus for each pair depending on */
+ /* compensatory/consistent mutations and incompatible seqs */
+ /* should be 0 for conserved pairs, >0 for good pairs */
+#define NONE -10000 /* score for forbidden pairs */
+ int n,i,j,k,l,s,score;
+ int dm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
+ {0,0,2,2,1,2,2} /* CG */,
+ {0,2,0,1,2,2,2} /* GC */,
+ {0,2,1,0,2,1,2} /* GU */,
+ {0,1,2,2,0,2,1} /* UG */,
+ {0,2,2,1,2,0,2} /* AU */,
+ {0,2,2,2,1,2,0} /* UA */};
+ n=Sali[0][0]; /* length of seqs */
+ for (i=1; i<n; i++) {
+ for (j=i+1; (j<i+TURN+1) && (j<=n); j++)
+ pscore[indx[j]+i] = NONE;
+ for (j=i+TURN+1; j<=n; j++) {
+ int pfreq[8]={0,0,0,0,0,0,0,0};
+ for (s=0; s<n_seq; s++) {
+ int type;
+ if (Sali[s][i]==0 && Sali[s][j]==0) type = 7; /* gap-gap */
+ else {
+ if ((AS[s][i] == '~')||(AS[s][j] == '~')) type = 7;
+ else type = pair[Sali[s][i]][Sali[s][j]];
+ }
+
+ pfreq[type]++;
+ }
+ if (pfreq[0]*2>n_seq) { pscore[indx[j]+i] = NONE; continue;}
+ for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
+ for (l=k+1; l<=6; l++)
+ /* scores for replacements between pairtypes */
+ /* consistent or compensatory mutations score 1 or 2 */
+ score += pfreq[k]*pfreq[l]*dm[k][l];
+ /* counter examples score -1, gap-gap scores -0.25 */
+ pscore[indx[j]+i] = cv_fact *
+ ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
+ }
+ }
+
+ if (noLonelyPairs) /* remove unwanted pairs */
+ for (k=1; k<n-TURN-1; k++)
+ for (l=1; l<=2; l++) {
+ int type,ntype=0,otype=0;
+ i=k; j = i+TURN+l;
+ type = pscore[indx[j]+i];
+ while ((i>=1)&&(j<=n)) {
+ if ((i>1)&&(j<n)) ntype = pscore[indx[j+1]+i-1];
+ if ((otype<-4*UNIT)&&(ntype<-4*UNIT)) /* worse than 2 counterex */
+ pscore[indx[j]+i] = NONE; /* i.j can only form isolated pairs */
+ otype = type;
+ type = ntype;
+ i--; j++;
+ }
+ }
+
+
+ if (fold_constrained&&(structure!=NULL)) {
+ int psij, hx, hx2, *stack, *stack2;
+ stack = (int *) space(sizeof(int)*(n+1));
+ stack2 = (int *) space(sizeof(int)*(n+1));
+
+ for(hx=hx2=0, j=1; j<=n; j++) {
+ switch (structure[j-1]) {
+ case 'x': /* can't pair */
+ for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
+ for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
+ break;
+ case '(':
+ stack[hx++]=j;
+ /* fallthrough */
+ case '[':
+ stack2[hx2++]=j;
+ /* fallthrough */
+ case '<': /* pairs upstream */
+ for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
+ break;
+ case ']':
+ if (hx2<=0) {
+ fprintf(stderr, "%s\n", structure);
+ nrerror("unbalanced brackets in constraints");
+ }
+ i = stack2[--hx2];
+ pscore[indx[j]+i]=NONE;
+ break;
+ case ')':
+ if (hx<=0) {
+ fprintf(stderr, "%s\n", structure);
+ nrerror("unbalanced brackets in constraints");
+ }
+ i = stack[--hx];
+ psij = pscore[indx[j]+i]; /* store for later */
+ for (k=j; k<=n; k++)
+ for (l=i; l<=j; l++)
+ pscore[indx[k]+l] = NONE;
+ for (l=i; l<=j; l++)
+ for (k=1; k<=i; k++)
+ pscore[indx[l]+k] = NONE;
+ for (k=i+1; k<j; k++)
+ pscore[indx[k]+i] = pscore[indx[j]+k] = NONE;
+ pscore[indx[j]+i] = (psij>0) ? psij : 0;
+ /* fallthrough */
+ case '>': /* pairs downstream */
+ for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
+ break;
+ }
+ }
+ if (hx!=0) {
+ fprintf(stderr, "%s\n", structure);
+ nrerror("unbalanced brackets in constraint string");
+ }
+ free(stack); free(stack2);
+ }
+}
+
+float alisnofold(const char **strings, const int max_assym, const int threshloop,
+ const int min_s2, const int max_s2, const int half_stem, const int max_half_stem) {
+ int s,n_seq, length, energy;
+ char * structure;
+ length = (int) strlen(strings[0]);
+ /* structure = (char *) space((unsigned) length+1); */
+ structure = NULL;
+ if (length>init_length) alisnoinitialize_fold(length);
+ if (fabs(P->temperature - temperature)>1e-6) snoupdate_fold_params();
+ for (s=0; strings[s]!=NULL; s++);
+ n_seq = s;
+ Sali = (short **) space(n_seq*sizeof(short *));
+ for (s=0; s<n_seq; s++) {
+ if (strlen(strings[s]) != length) nrerror("uneqal seqence lengths");
+ Sali[s] = aliencode_seq(strings[s]);
+ }
+ make_pscores((const short **) Sali, (const char *const *) strings, n_seq, structure);
+ energy=alifill_arrays(strings, max_assym, threshloop, min_s2, max_s2, half_stem, max_half_stem);
+ alibacktrack((const char **)strings, 0);
+ structure = (char *) space((length+1)*sizeof(char));
+ parenthesis_structure(structure, length);
+
+ free(structure);
+ for (s=0; s<n_seq; s++) free(Sali[s]);
+ free(Sali);
+ /* free(structure); */
+ /* free(S)*/; free(S1); /* free(BP); */
+ return (float) energy/100.;
+}
+
+PRIVATE int alifill_arrays(const char **strings, const int max_asymm, const int threshloop,
+ const int min_s2, const int max_s2, const int half_stem,
+ const int max_half_stem) {
+
+ int i, j, length, energy;
+ /* int decomp, new_fML; */
+ int *type, type_2;
+ int bonus,n_seq,s;
+
+
+ for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
+ type = (int *) space(n_seq*sizeof(int));
+ length = strlen(strings[0]);
+ bonus=0;
+ /* max_separation = (int) ((1.-LOCALITY)*(double)(length-2));*/ /* not in use */
+
+ /* for (i=(j>TURN?(j-TURN):1); i<j; i++) { */
+ /* } */
+ for (i = (length)-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
+ for (j = i+TURN+1; j <= length; j++) {
+ int p, q, ij,psc;
+ ij = indx[j]+i;
+ for (s=0; s<n_seq; s++) {
+ type[s] = pair[Sali[s][i]][Sali[s][j]];
+ if (type[s]==0) type[s]=7;
+ }
+ psc = pscore[indx[j]+i];
+ if (psc>=MINPSCORE) { /* we have a pair */
+ int new_c=0, stackEnergy=INF; /* seems that new_c immer den minimum von cij enthaelt */
+ /* hairpin ----------------------------------------------*/
+
+ for (new_c=s=0; s<n_seq; s++)
+ new_c += E_Hairpin(j-i-1,type[s],Sali[s][i+1],Sali[s][j-1],strings[s]+i-1,P);
+ /*--------------------------------------------------------
+ check for elementary structures involving more than one
+ closing pair (interior loop).
+ --------------------------------------------------------*/
+
+ for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
+ int minq = j-i+p-MAXLOOP-2;
+ if (minq<p+1+TURN) minq = p+1+TURN;
+ for (q = minq; q < j; q++) {
+ if (pscore[indx[q]+p]<MINPSCORE) continue;
+ if(abs((p-i) - (j-q)) > max_asymm) continue;
+ for (energy = s=0; s<n_seq; s++) {
+ type_2 = pair[Sali[s][q]][Sali[s][p]]; /* q,p not p,q! */
+ if (type_2 == 0) type_2 = 7;
+ energy += E_IntLoop(p-i-1, j-q-1, type[s], type_2,
+ Sali[s][i+1], Sali[s][j-1],
+ Sali[s][p-1], Sali[s][q+1],P);
+ }
+ new_c = MIN2(energy+c[indx[q]+p], new_c);
+ if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
+
+ } /* end q-loop */
+ } /* end p-loop */
+
+ /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
+
+ new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
+ cc[j] = new_c - psc; /* add covariance bonnus/penalty */
+ c[ij]=cc[j];
+ } /* end >> if (pair) << */
+ else c[ij] = INF;
+ /* done with c[i,j], now compute fML[i,j] */
+ /* free ends ? -----------------------------------------*/
+
+ }
+
+ {
+ int *FF; /* rotate the auxilliary arrays */
+ FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
+ FF = cc1; cc1=cc; cc=FF;
+ for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
+ }
+ }
+ foldlist = (folden**) space((length)*sizeof(folden*));
+
+ for(i=0; i< length; i++){
+ foldlist[i]=(folden*) space(sizeof(folden));
+ foldlist[i]->next=NULL;
+ foldlist[i]->k=INF+1;
+ foldlist[i]->energy=INF;
+
+ }
+ folden* head; /* we save the stem loop information in a list like structure */
+
+ for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
+ int max_k, min_k;
+ max_k = MIN2(length-min_s2,i+max_half_stem+1);
+ min_k = MAX2(i+half_stem+1, length-max_s2);
+ for (j = i+TURN+1; j <= length; j++) {
+ int ij,a,b;
+ ij = indx[j]+i;
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ mLoop[ij]=MIN2(mLoop[ij], c[indx[j-a]+i+b]);
+
+ }
+ }
+ if(mLoop[ij]>=n_seq*threshloop){
+ mLoop[ij]=INF;
+ }
+ else{
+ if(j>=min_k-1 && j < max_k){ /* comment if out to recover the known behaviour */
+ head = (folden*) space(sizeof(folden));
+ head->k=j;
+ head->energy=mLoop[ij];
+ head->next=foldlist[i];
+ foldlist[i] = head;
+
+ }
+ }
+ }
+
+ }
+ free(type);
+ return mLoop[indx[length]+1];/* mLoop; */
+}
+
+PRIVATE int alibacktrack(const char **strings, int s) {
+
+ /*------------------------------------------------------------------
+ trace back through the "c", "f5" and "fML" arrays to get the
+ base pairing list. No search for equivalent structures is done.
+ This is fast, since only few structure elements are recalculated.
+ ------------------------------------------------------------------*/
+
+ /* normally s=0.
+ If s>0 then s items have been already pushed onto the sector stack */
+ int i, j, length, energy;/* , new; */
+ int type_2;
+ int bonus,n_seq,*type; int b=0,cov_en = 0;
+
+ length = strlen(strings[0]);
+ for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
+ type = (int *) space(n_seq*sizeof(int));
+ if (s==0) {
+ sector[++s].i = 1;
+ sector[s].j = length;
+ sector[s].ml = 2 ;
+ }
+ while (s>0) {
+ int ml, ss, cij, traced, i1, j1, p, q;
+ int canonical = 1; /* (i,j) closes a canonical structure */
+ i = sector[s].i;
+ j = sector[s].j;
+ ml = sector[s--].ml; /* ml is a flag indicating if backtracking is to
+ occur in the fML- (1) or in the f-array (0) */
+ if (ml==2) {
+ base_pair[++b].i = i;
+ base_pair[b].j = j;
+ goto repeat1;
+ }
+
+ if (j < i+TURN+1) continue; /* no more pairs in this interval */
+
+
+ repeat1:
+
+ /*----- begin of "repeat:" -----*/
+ if (canonical) cij = c[indx[j]+i];
+ for (ss=0; ss<n_seq; ss++) {
+ type[ss] = pair[Sali[ss][i]][Sali[ss][j]];
+ if (type[ss]==0) type[ss] = 7;
+ }
+ bonus = 0;
+
+ if (noLonelyPairs)
+ if (cij == c[indx[j]+i]) {
+ /* (i.j) closes canonical structures, thus
+ (i+1.j-1) must be a pair */
+ for (ss=0; ss<n_seq; ss++) {
+ type_2 = pair[Sali[ss][j-1]][Sali[ss][i+1]]; /* j,i not i,j */
+ if (type_2==0) type_2 = 7;
+ cij -= P->stack[type[ss]][type_2];
+ }
+ cij += pscore[indx[j]+i];
+ base_pair[++b].i = i+1;
+ base_pair[b].j = j-1;
+ cov_en += pscore[indx[j-1]+i+1];
+ i++; j--;
+ canonical=0;
+ goto repeat1;
+ }
+ canonical = 1;
+ cij += pscore[indx[j]+i];
+ {int cc=0;
+ for (ss=0; ss<n_seq; ss++)
+ cc += E_Hairpin(j-i-1, type[ss], Sali[ss][i+1], Sali[ss][j-1], strings[ss]+i-1,P);
+ if (cij == cc) /* found hairpin */
+ continue;
+ }
+ for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
+ int minq;
+ minq = j-i+p-MAXLOOP-2;
+ if (minq<p+1+TURN) minq = p+1+TURN;
+ for (q = j-1; q >= minq; q--) {
+ for (ss=energy=0; ss<n_seq; ss++) {
+ type_2 = pair[Sali[ss][q]][Sali[ss][p]]; /* q,p not p,q */
+ if (type_2==0) type_2 = 7;
+ energy += E_IntLoop(p-i-1, j-q-1, type[ss], type_2,
+ Sali[ss][i+1], Sali[ss][j-1],
+ Sali[ss][p-1], Sali[ss][q+1],P);
+ }
+ traced = (cij == energy+c[indx[q]+p]);
+ if (traced) {
+ base_pair[++b].i = p;
+ base_pair[b].j = q;
+ cov_en += pscore[indx[q]+p];
+ i = p, j = q;
+ goto repeat1;
+ }
+ }
+ }
+
+ /* end of repeat: --------------------------------------------------*/
+
+ /* (i.j) must close a multi-loop */
+ /* tt = rtype[type]; */
+/* mm = bonus+P->MLclosing+P->MLintern[tt]; */
+/* d5 = P->dangle5[tt][S1[j-1]]; */
+/* d3 = P->dangle3[tt][S1[i+1]]; */
+ i1 = i+1; j1 = j-1;
+ sector[s+1].ml = sector[s+2].ml = 1;
+
+/* if (k<=j-3-TURN) { /\\* found the decomposition *\\/ *\/ */
+/* sector[++s].i = i1; */
+/* sector[s].j = k; */
+/* sector[++s].i = k+1; */
+/* sector[s].j = j1; */
+/* } /\* else { *\/ */
+/* nrerror("backtracking failed in repeat"); */
+/* } */
+
+ }
+ base_pair[0].i = b; /* save the total number of base pairs */
+ free(type);
+ return cov_en;
+}
+
+PRIVATE int fill_arrays(const char *string, const int max_asymm, const int threshloop,
+ const int min_s2, const int max_s2, const int half_stem, const int max_half_stem) {
+
+ int i, j, length, energy;
+ /* int decomp;*/ /*, new_fML; */
+ int no_close, type, type_2;
+ int bonus;
+ int min_c;
+
+ min_c=INF;
+ length = (int) strlen(string);
+ bonus=0;
+ /* max_separation = (int) ((1.-LOCALITY)*(double)(length-2)); */ /* not in use */
+
+
+
+
+ for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
+ /* printf("i=%d\t",i); */
+ for (j = i+TURN+1; j <= length; j++) {
+/* printf("j=%d,",j); */
+ int p, q, ij;
+ ij = indx[j]+i;
+ type = ptype[ij];
+ bonus = 0;
+ energy = INF;
+
+ if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
+ if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
+ if ((BP[i]==-4)||(BP[j]==-4)) type=0;
+
+ no_close = (((type==3)||(type==4))&&no_closingGU);
+
+ /* if (j-i-1 > max_separation) type = 0; */ /* forces locality degree */
+
+ if (type) { /* we have a pair */
+ int new_c=0, stackEnergy=INF; /* seems that new_c immer den minimum von cij enthaelt */
+ /* hairpin ----------------------------------------------*/
+
+ if (no_close) new_c = FORBIDDEN;
+ else
+ new_c = E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], string+i-1,P); /* computes hair pin structure for subsequence i...j */
+
+ /*--------------------------------------------------------
+ check for elementary structures involving more than one
+ closing pair (interior loop).
+ --------------------------------------------------------*/
+
+ for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
+ int minq = j-i+p-MAXLOOP-2;
+ if (minq<p+1+TURN) minq = p+1+TURN;
+ for (q = minq; q < j; q++) {
+
+ if(abs((p-i) - (j-q)) > max_asymm) continue;
+ type_2 = ptype[indx[q]+p];
+
+ if (type_2==0) continue;
+ type_2 = rtype[type_2];
+
+ if (no_closingGU)
+ if (no_close||(type_2==3)||(type_2==4))
+ if ((p>i+1)||(q<j-1)) continue; /* continue unless stack */
+
+ energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
+ S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
+ new_c = MIN2(energy+c[indx[q]+p], new_c);
+ if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
+
+ } /* end q-loop */
+ } /* end p-loop */
+
+
+
+
+ /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
+
+
+ new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
+ cc[j] = new_c;
+ c[ij] = new_c;
+ /* min_c=MIN2(min_c, c[ij]); */
+
+ } /* end >> if (pair) << */
+
+ else c[ij] = INF;
+
+
+ /* done with c[i,j], now compute fML[i,j] */
+ /* free ends ? -----------------------------------------*/
+
+ }
+
+ {
+ int *FF; /* rotate the auxilliary arrays */
+ FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
+ FF = cc1; cc1=cc; cc=FF;
+ for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
+ }
+ }
+ foldlist = (folden**) space((length+1)*sizeof(folden*));
+ foldlist_XS = (folden**) space((length+1)*sizeof(folden*));
+ /* linked list initialization*/
+ for(i=0; i<=length; i++){
+ foldlist[i]=(folden*) space(sizeof(folden));
+ foldlist[i]->next=NULL;
+ foldlist[i]->k=INF+1;
+ foldlist[i]->energy=INF;
+ foldlist_XS[i]=(folden*) space(sizeof(folden));
+ foldlist_XS[i]->next=NULL;
+ foldlist_XS[i]->k=INF+1;
+ foldlist_XS[i]->energy=INF;
+ }
+ folden* head; /* we save the stem loop information in a list like structure */
+ folden* head_XS;
+ for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
+ int max_k, min_k;
+ max_k = MIN2(length-min_s2,i+max_half_stem+1);
+ min_k = MAX2(i+half_stem+1, length-max_s2);
+
+
+ for (j = i+TURN+1; j <= length; j++) {
+ int ij,a,b;
+ ij = indx[j]+i;
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ mLoop[ij]=MIN2(mLoop[ij], c[indx[j-a]+i+b]);
+ /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i]); */
+ /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j]+i+1]); */
+ /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-1]+i+1]); */
+ /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i+1]); */
+ /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j]+i+2]); */
+ /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-1]+i+2]); */
+ /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i+2]); */
+ }
+ }
+ min_c = MIN2(mLoop[ij] ,min_c);
+
+ if(mLoop[ij]>=threshloop){
+ mLoop[ij]=INF;
+ }
+ else{
+ if(j>=min_k-1 && j <= max_k){ /* comment if out to recover the known behaviour */
+ head = (folden*) space(sizeof(folden));
+ head->k=j;
+ head->energy=mLoop[ij];
+ head->next=foldlist[i];
+ foldlist[i] = head;
+ head_XS = (folden*) space(sizeof(folden));
+ head_XS->k=i;
+ head_XS->energy=mLoop[ij];
+ head_XS->next=foldlist_XS[j];
+ foldlist_XS[j] = head_XS;
+ }
+ }
+ }
+
+ }
+/* int count=0; */
+/* for(i=0; i< length; i++){ */
+/* folden *temp; */
+/* temp = foldlist[i]; */
+/* while(temp->next){ */
+/* count++; */
+/* printf("count %d: i%d j%d energy %d \n", count, i, temp->k, temp->energy); */
+/* temp=temp->next; */
+/* } */
+/* } */
+/* printf("Count %d \n", count); */
+/* count=0; */
+/* for(i=length-1; i>=0; i--){ */
+/* folden *temp; */
+/* temp = foldlist_XS[i]; */
+/* while(temp->next){ */
+/* count++; */
+/* printf("count %d: i%d j%d energy %d \n", count, temp->k,i, temp->energy); */
+/* temp=temp->next; */
+/* } */
+/* } */
+/* printf("Count %d \n", count); */
+/* return mLoop[indx[length]+1]; */ /* mLoop; */
+ return min_c;
+ /* printf("\nmin_array = %d\n", min_c); */
+ /* return f5[length]; */
+}
+
+
+
+
+PRIVATE void backtrack(const char *string, int s) {
+
+ /*------------------------------------------------------------------
+ trace back through the "c", "f5" and "fML" arrays to get the
+ base pairing list. No search for equivalent structures is done.
+ This is fast, since only few structure elements are recalculated.
+ ------------------------------------------------------------------*/
+
+ /* normally s=0.
+ If s>0 then s items have been already pushed onto the sector stack */
+ int i, j, /*k,*/ length, energy, new;
+ int no_close, type, type_2;/* , tt; */
+ int bonus;
+ int b=0;
+
+ length = strlen(string);
+ if (s==0) {
+ sector[++s].i = 1;
+ sector[s].j = length;
+ sector[s].ml = 2 ;
+ }
+ while (s>0) {
+ int ml, cij, traced, i1, j1, /*d3, d5, mm,*/ p, q;
+ int canonical = 1; /* (i,j) closes a canonical structure */
+ i = sector[s].i;
+ j = sector[s].j;
+ ml = sector[s--].ml; /* ml is a flag indicating if backtracking is to
+ occur in the fML- (1) or in the f-array (0) */
+ if (ml==2) {
+ base_pair[++b].i = i;
+ base_pair[b].j = j;
+ goto repeat1;
+ }
+
+ if (j < i+TURN+1) continue; /* no more pairs in this interval */
+
+
+ repeat1:
+
+ /*----- begin of "repeat:" -----*/
+ if (canonical) cij = c[indx[j]+i];
+ type = ptype[indx[j]+i];
+ bonus = 0;
+ if (fold_constrained) {
+ if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
+ if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
+ }
+ if (noLonelyPairs)
+ if (cij == c[indx[j]+i]) {
+ /* (i.j) closes canonical structures, thus
+ (i+1.j-1) must be a pair */
+ type_2 = ptype[indx[j-1]+i+1]; type_2 = rtype[type_2];
+ cij -= P->stack[type][type_2] + bonus;
+ base_pair[++b].i = i+1;
+ base_pair[b].j = j-1;
+ i++; j--;
+ canonical=0;
+ goto repeat1;
+ }
+ canonical = 1;
+ no_close = (((type==3)||(type==4))&&no_closingGU&&(bonus==0));
+ if (no_close) {
+ if (cij == FORBIDDEN) continue;
+ } else
+ if (cij == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1],string+i-1,P)+bonus)
+ continue;
+ for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
+ int minq;
+ minq = j-i+p-MAXLOOP-2;
+ if (minq<p+1+TURN) minq = p+1+TURN;
+ for (q = j-1; q >= minq; q--) {
+ type_2 = ptype[indx[q]+p];
+ if (type_2==0) continue;
+ type_2 = rtype[type_2];
+ if (no_closingGU)
+ if (no_close||(type_2==3)||(type_2==4))
+ if ((p>i+1)||(q<j-1)) continue; /* continue unless stack */
+ energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
+ S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
+ new = energy+c[indx[q]+p]+bonus;
+ traced = (cij == new);
+ if (traced) {
+ base_pair[++b].i = p;
+ base_pair[b].j = q;
+ i = p, j = q;
+ goto repeat1;
+ }
+ }
+ }
+
+ /* end of repeat: --------------------------------------------------*/
+
+ /* (i.j) must close a multi-loop */
+/* tt = rtype[type]; */
+/* mm = bonus+P->MLclosing+P->MLintern[tt]; */
+/* d5 = P->dangle5[tt][S1[j-1]]; */
+/* d3 = P->dangle3[tt][S1[i+1]]; */
+ i1 = i+1; j1 = j-1;
+ sector[s+1].ml = sector[s+2].ml = 1;
+
+/* if (k<=j-3-TURN) { */ /* found the decomposition */
+/* sector[++s].i = i1; */
+/* sector[s].j = k; */
+/* sector[++s].i = k+1; */
+/* sector[s].j = j1; */
+/* } else { */
+/* nrerror("backtracking failed in repeat"); */
+/* } */
+/* */
+ }
+
+ base_pair[0].i = b; /* save the total number of base pairs */
+}
+
+char *snobacktrack_fold_from_pair(const char *sequence, int i, int j) {
+ char *structure;
+ sector[1].i = i;
+ sector[1].j = j;
+ sector[1].ml = 2;
+ base_pair[0].i=0;
+ encode_seq(sequence);
+ backtrack(sequence, 1);
+ structure = (char *) space((strlen(sequence)+1)*sizeof(char));
+ parenthesis_structure(structure, strlen(sequence));
+ free(S);free(S1);
+ return structure;
+}
+
+char *alisnobacktrack_fold_from_pair(const char **strings, int i, int j, int *cov) {
+ char *structure;
+ int n_seq, s, length;
+ length = (int) strlen(strings[0]);
+ for (s=0; strings[s]!=NULL; s++);
+ n_seq = s;
+ sector[1].i = i;
+ sector[1].j = j;
+ sector[1].ml = 2;
+ base_pair[0].i=0;
+ /* encode_seq(sequence); */
+ Sali = (short **) space(n_seq*sizeof(short *));
+ for (s=0; s<n_seq; s++) {
+ if (strlen(strings[s]) != length) nrerror("uneqal seqence lengths");
+ Sali[s] = aliencode_seq(strings[s]);
+ }
+ *cov=alibacktrack(strings, 1);
+ structure = (char *) space((length+1)*sizeof(char));
+ parenthesis_structure(structure, length);
+ free(S);free(S1);
+ for (s=0; s<n_seq; s++) {
+ free(Sali[s]);
+ }
+ free(Sali);
+ return structure;
+}
+
+
+
+/*---------------------------------------------------------------------------*/
+
+
+/*---------------------------------------------------------------------------*/
+
+
+/*--------------------------------------------------------------------------*/
+
+
+/*---------------------------------------------------------------------------*/
+
+
+PRIVATE void encode_seq(const char *sequence) {
+ unsigned int i,l;
+
+ l = strlen(sequence);
+ S = (short *) space(sizeof(short)*(l+2));
+ S1= (short *) space(sizeof(short)*(l+2));
+ /* S1 exists only for the special X K and I bases and energy_set!=0 */
+ S[0] = (short) l;
+
+ for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
+ S[i]= (short) encode_char(toupper(sequence[i-1]));
+ S1[i] = alias[S[i]]; /* for mismatches of nostandard bases */
+ }
+ /* for circular folding add first base at position n+1 and last base at
+ position 0 in S1 */
+ S[l+1] = S[1]; S1[l+1]=S1[1]; S1[0] = S1[l];
+}
+
+PRIVATE short * aliencode_seq(const char *sequence) {
+ unsigned int i,l;
+ short *Stemp;
+ l = strlen(sequence);
+ Stemp = (short *) space(sizeof(short)*(l+2));
+ Stemp[0] = (short) l;
+
+ /* make numerical encoding of sequence */
+ for (i=1; i<=l; i++)
+ Stemp[i]= (short) encode_char(toupper(sequence[i-1]));
+
+ /* for circular folding add first base at position n+1 */
+ /* Stemp[l+1] = Stemp[1]; */
+
+ return Stemp;
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void letter_structure(char *structure, int length)
+{
+ int n, k, x, y;
+
+ for (n = 0; n <= length-1; structure[n++] = ' ') ;
+ structure[length] = '\0';
+
+ for (n = 0, k = 1; k <= base_pair[0].i; k++) {
+ y = base_pair[k].j;
+ x = base_pair[k].i;
+ if (x-1 > 0 && y+1 <= length) {
+ if (structure[x-2] != ' ' && structure[y] == structure[x-2]) {
+ structure[x-1] = structure[x-2];
+ structure[y-1] = structure[x-1];
+ continue;
+ }
+ }
+ if (structure[x] != ' ' && structure[y-2] == structure[x]) {
+ structure[x-1] = structure[x];
+ structure[y-1] = structure[x-1];
+ continue;
+ }
+ n++;
+ structure[x-1] = alpha[n-1];
+ structure[y-1] = alpha[n-1];
+ }
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void parenthesis_structure(char *structure, int length)
+{
+ int n, k;
+
+ for (n = 0; n <= length-1; structure[n++] = '.') ;
+ structure[length] = '\0';
+
+ for (k = 1; k <= base_pair[0].i; k++) {
+ structure[base_pair[k].i-1] = '(';
+ structure[base_pair[k].j-1] = ')';
+ }
+}
+/*---------------------------------------------------------------------------*/
+
+PUBLIC void snoupdate_fold_params(void)
+{
+ if(P) free(P);
+ P = scale_parameters();
+ make_pair_matrix();
+ if (init_length < 0) init_length=0;
+}
+
+/*---------------------------------------------------------------------------*/
+/* PRIVATE short *pair_table; */
+
+
+/*---------------------------------------------------------------------------*/
+#if 0
+PRIVATE int stack_energy(int i, const char *string)
+{
+ /* calculate energy of substructure enclosed by (i,j) */
+ int ee, energy = 0;
+ int j, p, q, type;
+
+ j=pair_table[i];
+ type = pair[S[i]][S[j]];
+ if (type==0) {
+ type=7;
+ }
+
+ p=i; q=j;
+ while (p<q) { /* process all stacks and interior loops */
+ int type_2;
+ while (pair_table[++p]==0);
+ while (pair_table[--q]==0);
+ if ((pair_table[q]!=(short)p)||(p>q)) break;
+ type_2 = pair[S[q]][S[p]];
+ if (type_2==0) {
+ type_2=7;
+ }
+ /* energy += LoopEnergy(i, j, p, q, type, type_2); */
+ if ( SAME_STRAND(i,p) && SAME_STRAND(q,j) )
+ ee = LoopEnergy(p-i-1, j-q-1, type, type_2,
+ S1[i+1], S1[j-1], S1[p-1], S1[q+1]);
+ energy += ee;
+ i=p; j=q; type = rtype[type_2];
+ } /* end while */
+
+ /* p,q don't pair must have found hairpin or multiloop */
+
+ if (p>q) { /* hair pin */
+ if (SAME_STRAND(i,j))
+ ee = snoHairpinE(j-i-1, type, S1[i+1], S1[j-1], string+i-1);
+ energy += ee;
+
+ return energy;
+ }
+
+ /* (i,j) is exterior pair of multiloop */
+ while (p<j) {
+ /* add up the contributions of the substructures of the ML */
+ energy += stack_energy(p, string);
+ p = pair_table[p];
+ /* search for next base pair in multiloop */
+ while (pair_table[++p]==0);
+ }
+ {
+ int ii;
+ ii = cut_in_loop(i);
+ }
+ energy += ee;
+
+ return energy;
+}
+
+/*---------------------------------------------------------------------------*/
+
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE int cut_in_loop(int i) {
+ /* walk around the loop; return j pos of pair after cut if
+ cut_point in loop else 0 */
+ int p, j;
+ p = j = pair_table[i];
+ do {
+ i = pair_table[p]; p = i+1;
+ while ( pair_table[p]==0 ) p++;
+ } while (p!=j && SAME_STRAND(i,p));
+ return SAME_STRAND(i,p) ? 0 : pair_table[p];
+}
+#endif
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void make_ptypes(const short *S, const char *structure) {
+ int n,i,j,k,l;
+
+ n=S[0];
+ for (k=1; k<n-TURN; k++)
+ for (l=1; l<=2; l++) {
+ int type,ntype=0,otype=0;
+ i=k; j = i+TURN+l; 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[indx[j]+i] = (char) type;
+ otype = type;
+ type = ntype;
+ i--; j++;
+ }
+ }
+
+ if (fold_constrained&&(structure!=NULL)) {
+ constrain_ptypes(structure, (unsigned int)n, ptype, BP, TURN, 0);
+#if 0
+ int hx, *stack;
+ char type;
+ stack = (int *) space(sizeof(int)*(n+1));
+
+ for(hx=0, j=1; j<=n; j++) {
+ switch (structure[j-1]) {
+ case '|': BP[j] = -1; break;
+ case 'x': /* can't pair */
+ for (l=1; l<j-TURN; l++) ptype[indx[j]+l] = 0;
+ for (l=j+TURN+1; l<=n; l++) ptype[indx[l]+j] = 0;
+ break;
+ case '(':
+ stack[hx++]=j;
+ /* fallthrough */
+ case '<': /* pairs upstream */
+ for (l=1; l<j-TURN; l++) ptype[indx[j]+l] = 0;
+ break;
+ case ')':
+ if (hx<=0) {
+ fprintf(stderr, "%s\n", structure);
+ nrerror("unbalanced brackets in constraints");
+ }
+ i = stack[--hx];
+ type = ptype[indx[j]+i];
+ for (k=i+1; k<=n; k++) ptype[indx[k]+i] = 0;
+ /* don't allow pairs i<k<j<l */
+ for (l=j; l<=n; l++)
+ for (k=i+1; k<=j; k++) ptype[indx[l]+k] = 0;
+ /* don't allow pairs k<i<l<j */
+ for (l=i; l<=j; l++)
+ for (k=1; k<=i; k++) ptype[indx[l]+k] = 0;
+ for (k=1; k<j; k++) ptype[indx[j]+k] = 0;
+ ptype[indx[j]+i] = (type==0)?7:type;
+ /* fallthrough */
+ case '>': /* pairs downstream */
+ for (l=j+TURN+1; l<=n; l++) ptype[indx[l]+j] = 0;
+ break;
+ }
+ }
+ if (hx!=0) {
+ fprintf(stderr, "%s\n", structure);
+ nrerror("unbalanced brackets in constraint string");
+ }
+ free(stack);
+#endif
+ }
+}
git://source.jalview.org / jabaws.git / blobdiff