4 c Ivo L Hofacker and Walter Fontana
7 /* Last changed Time-stamp: <2008-11-25 16:34:36 ivo> */
15 #include "../config.h"
22 #define PRIVATE static
26 PUBLIC unsigned short xsubi[3];
28 PRIVATE char scale1[] = "....,....1....,....2....,....3....,....4";
29 PRIVATE char scale2[] = "....,....5....,....6....,....7....,....8";
30 PRIVATE char *inbuf = NULL;
32 PRIVATE char *inbuf2 = NULL;
33 PRIVATE unsigned int typebuf2 = 0;
35 /*-------------------------------------------------------------------------*/
37 PUBLIC void *space(unsigned size) {
40 if ( (pointer = (void *) calloc(1, (size_t) size)) == NULL) {
43 fprintf(stderr,"SPACE: requested size: %d\n", size);
44 nrerror("SPACE allocation failure -> EINVAL");
48 nrerror("SPACE allocation failure -> no memory");
54 #define space(S) calloc(1,(S))
58 /* dmalloc.h #define's xrealloc */
59 void *xrealloc (void *p, unsigned size) {
62 p = (void *) realloc(p, size);
66 fprintf(stderr,"xrealloc: requested size: %d\n", size);
67 nrerror("xrealloc allocation failure -> EINVAL");
71 nrerror("xrealloc allocation failure -> no memory");
76 /*------------------------------------------------------------------------*/
78 PUBLIC void nrerror(const char message[]) /* output message upon error */
80 fprintf(stderr, "ERROR: %s\n", message);
84 PUBLIC void warn_user(const char message[]){
85 fprintf(stderr, "WARNING: %s\n", message);
88 /*------------------------------------------------------------------------*/
89 PUBLIC void init_rand(void)
93 xsubi[0] = xsubi[1] = xsubi[2] = (unsigned short) t; /* lower 16 bit */
94 xsubi[1] += (unsigned short) ((unsigned)t >> 6);
95 xsubi[2] += (unsigned short) ((unsigned)t >> 12);
97 srand((unsigned int) t);
101 /*------------------------------------------------------------------------*/
103 PUBLIC double urn(void)
104 /* uniform random number generator; urn() is in [0,1] */
105 /* uses a linear congruential library routine */
106 /* 48 bit arithmetic */
109 extern double erand48(unsigned short[]);
110 return erand48(xsubi);
112 return ((double) rand())/RAND_MAX;
116 /*------------------------------------------------------------------------*/
118 PUBLIC int int_urn(int from, int to)
120 return ( ( (int) (urn()*(to-from+1)) ) + from );
123 /*------------------------------------------------------------------------*/
125 PUBLIC void filecopy(FILE *from, FILE *to)
129 while ((c = getc(from)) != EOF) (void)putc(c, to);
132 /*-----------------------------------------------------------------*/
134 PUBLIC char *time_stamp(void)
138 cal_time = time(NULL);
139 return ( ctime(&cal_time) );
142 /*-----------------------------------------------------------------*/
144 PUBLIC char *random_string(int l, const char symbols[])
149 base = (int) strlen(symbols);
150 r = (char *) space(sizeof(char)*(l+1));
152 for (i = 0; i < l; i++) {
153 rn = (int) (urn()*base); /* [0, base-1] */
160 /*-----------------------------------------------------------------*/
162 PUBLIC int hamming(const char *s1, const char *s2)
166 for (; *s1 && *s2; s1++, s2++)
171 PUBLIC int hamming_bound(const char *s1, const char *s2, int boundary)
175 for (; *s1 && *s2 && boundary; s1++, s2++, boundary--)
179 /*-----------------------------------------------------------------*/
181 PUBLIC char *get_line(FILE *fp) /* reads lines of arbitrary length from fp */
183 char s[512], *line, *cp;
184 int len=0, size=0, l;
187 if (fgets(s, 512, fp)==NULL) break;
188 cp = strchr(s, '\n');
189 if (cp != NULL) *cp = '\0';
190 l = len + (int)strlen(s);
192 size = (int)((l+1)*1.2);
193 line = (char *) xrealloc(line, size*sizeof(char));
202 PUBLIC int skip_comment_lines(char **line){
203 if((*line = get_line(stdin))==NULL) return -1;
205 while((**line=='*')||(**line=='\0')){
207 if((*line = get_line(stdin))==NULL) return -1;
212 PUBLIC unsigned int get_input_line(char **string, unsigned int option){
217 * read lines until informative data appears or
218 * report an error if anything goes wrong
220 if((line = get_line(stdin))==NULL) return VRNA_INPUT_ERROR;
222 if(!(option & VRNA_INPUT_NOSKIP_COMMENTS))
223 while((*line=='*')||(*line=='\0')){
225 if((line = get_line(stdin))==NULL) return VRNA_INPUT_ERROR;
228 l = (int) strlen(line);
230 /* break on '@' sign if not disabled */
233 return VRNA_INPUT_QUIT;
235 /* print line read if not disabled */
236 /* if(!(option & VRNA_INPUT_NOPRINT)) printf("%s\n", line); */
238 /* eliminate whitespaces at the end of the line read */
239 if(!(option & VRNA_INPUT_NO_TRUNCATION)){
240 for(i = l-1; i >= 0; i--){
241 if (line[i] == ' ') continue;
242 else if (line[i] == '\t') continue;
245 line[(i >= 0) ? (i+1) : 0] = '\0';
250 /* alloc memory for the string */
251 *string = (char *) space(sizeof(char) * (strlen(line) + 1));
252 r = VRNA_INPUT_FASTA_HEADER;
253 i = sscanf(line, ">%s", *string);
255 i = (int) strlen(*string);
256 *string = (char *) xrealloc(*string, (i+1)*sizeof(char));
264 return VRNA_INPUT_ERROR;
268 *string = strdup(line);
271 return VRNA_INPUT_MISC;
274 PUBLIC unsigned int get_multi_input_line(char **string, unsigned int option){
280 line = (inbuf) ? inbuf : get_line(stdin);
285 * read lines until informative data appears or
286 * report an error if anything goes wrong
288 if(!line) return VRNA_INPUT_ERROR;
290 l = (int)strlen(line);
292 /* eliminate whitespaces at the end of the line read */
293 if(!(option & VRNA_INPUT_NO_TRUNCATION)){
294 for(i = l-1; i >= 0; i--){
295 if (line[i] == ' ') continue;
296 else if (line[i] == '\t') continue;
299 line[(i >= 0) ? (i+1) : 0] = '\0';
302 l = (int)strlen(line);
303 str_length = (*string) ? (int) strlen(*string) : 0;
306 case '@': /* user abort */
307 if(state) inbuf = line;
309 return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_QUIT;
311 case '\0': /* empty line */
312 if(option & VRNA_INPUT_NOSKIP_BLANK_LINES){
313 if(state) inbuf = line;
315 return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_BLANK_LINE;
319 case '#': case '%': case ';': case '/': case '*': case ' ':
321 if(option & VRNA_INPUT_NOSKIP_COMMENTS){
322 if(state) inbuf = line;
324 return (state == 2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_COMMENT;
328 case '>': /* fasta header */
329 if(state) inbuf = line;
331 return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_FASTA_HEADER;
333 case 'x': case 'e': case 'l': case '&': /* seems to be a constraint or line starting with second sequence for dimer calculations */
335 /* lets see if this assumption holds for the complete line */
336 while((line[i] == 'x') || (line[i] == 'e') || (line[i] == 'l')) i++;
337 /* lines solely consisting of 'x's, 'e's or 'l's will be considered as structure constraint */
340 ((line[i]>64) && (line[i]<91)) /* A-Z */
341 || ((line[i]>96) && (line[i]<123)) /* a-z */
343 if(option & VRNA_INPUT_FASTA_HEADER){
344 /* are we in structure mode? Then we remember this line for the next round */
345 if(state == 2){ inbuf = line; return VRNA_INPUT_CONSTRAINT;}
347 *string = (char *)xrealloc(*string, sizeof(char) * (str_length + l + 1));
348 strcpy(*string + str_length, line);
353 /* otherwise return line read */
354 else{ *string = line; return VRNA_INPUT_SEQUENCE;}
356 /* mmmh? it really seems to be a constraint */
358 case '<': case '.': case '|': case '(': case ')': case '[': case ']': case '{': case '}': case ',': case '+':
359 /* seems to be a structure or a constraint */
360 /* either we concatenate this line to one that we read previously */
361 if(option & VRNA_INPUT_FASTA_HEADER){
364 return VRNA_INPUT_SEQUENCE;
367 *string = (char *)xrealloc(*string, sizeof(char) * (str_length + l + 1));
368 strcpy(*string + str_length, line);
372 /* or we return it as it is */
375 return VRNA_INPUT_CONSTRAINT;
378 default: if(option & VRNA_INPUT_FASTA_HEADER){
379 /* are we already in sequence mode? */
382 return VRNA_INPUT_CONSTRAINT;
385 *string = (char *)xrealloc(*string, sizeof(char) * (str_length + l + 1));
386 strcpy(*string + str_length, line);
390 /* otherwise return line read */
393 return VRNA_INPUT_SEQUENCE;
397 line = get_line(stdin);
400 return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_ERROR;
403 PUBLIC unsigned int read_record(char **header, char **sequence, char ***rest, unsigned int options){
404 unsigned int input_type, return_type, tmp_type;
409 return_type = tmp_type = 0;
410 input_string = *header = *sequence = NULL;
411 *rest = (char **)space(sizeof(char *));
413 /* remove unnecessary option flags from options variable... */
414 options &= ~VRNA_INPUT_FASTA_HEADER;
416 /* read first input or last buffered input */
418 input_type = typebuf2;
419 input_string = inbuf2;
423 else input_type = get_multi_input_line(&input_string, options);
425 if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) return input_type;
427 /* skip everything until we read either a fasta header or a sequence */
428 while(input_type & (VRNA_INPUT_MISC | VRNA_INPUT_CONSTRAINT | VRNA_INPUT_BLANK_LINE)){
429 free(input_string); input_string = NULL;
430 input_type = get_multi_input_line(&input_string, options);
431 if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) return input_type;
434 if(input_type & VRNA_INPUT_FASTA_HEADER){
435 return_type |= VRNA_INPUT_FASTA_HEADER; /* remember that we've read a fasta header */
436 *header = input_string;
438 /* get next data-block with fasta support if not explicitely forbidden by VRNA_INPUT_NO_SPAN */
439 input_type = get_multi_input_line(
441 ((options & VRNA_INPUT_NO_SPAN) ? 0 : VRNA_INPUT_FASTA_HEADER) | options
443 if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) return (return_type | input_type);
446 if(input_type & VRNA_INPUT_SEQUENCE){
447 return_type |= VRNA_INPUT_SEQUENCE; /* remember that we've read a sequence */
448 *sequence = input_string;
450 } else nrerror("sequence input missing");
452 /* read the rest until we find user abort, EOF, new sequence or new fasta header */
453 if(!(options & VRNA_INPUT_NO_REST)){
454 options |= VRNA_INPUT_NOSKIP_COMMENTS; /* allow commetns to appear in rest output */
455 tmp_type = VRNA_INPUT_QUIT | VRNA_INPUT_ERROR | VRNA_INPUT_SEQUENCE | VRNA_INPUT_FASTA_HEADER;
456 if(options & VRNA_INPUT_NOSKIP_BLANK_LINES) tmp_type |= VRNA_INPUT_BLANK_LINE;
457 while(!((input_type = get_multi_input_line(&input_string, options)) & tmp_type)){
458 *rest = xrealloc(*rest, sizeof(char **)*(++rest_count + 1));
459 (*rest)[rest_count-1] = input_string;
463 if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) return input_type;
466 /* finished reading everything...
467 * we now put the last line into the buffer if necessary
468 * since it should belong to the next record
470 inbuf2 = input_string;
471 typebuf2 = input_type;
473 (*rest)[rest_count] = NULL;
474 return (return_type);
478 /*-----------------------------------------------------------------*/
480 PUBLIC char *pack_structure(const char *struc) {
481 /* 5:1 compression using base 3 encoding */
483 unsigned char *packed;
485 l = (int) strlen(struc);
486 packed = (unsigned char *) space(((l+4)/5+1)*sizeof(unsigned char));
491 for (p=pi=0; pi<5; pi++) {
503 default: nrerror("pack_structure: illegal charcter in structure");
507 packed[j++] = (unsigned char) (p+1); /* never use 0, so we can use
510 packed[j] = '\0'; /* for str*() functions */
511 return (char *) packed;
514 PUBLIC char *unpack_structure(const char *packed) {
515 /* 5:1 compression using base 3 encoding */
518 unsigned const char *pp;
519 char code[3] = {'(', '.', ')'};
521 l = (int) strlen(packed);
522 pp = (const unsigned char *) packed;
523 struc = (char *) space((l*5+1)*sizeof(char)); /* up to 4 byte extra */
525 for (i=j=0; i<l; i++) {
526 register int p, c, k;
529 for (k=4; k>=0; k--) {
532 struc[j+k] = code[c];
537 while (struc[j] == '(') /* strip trailing ( */
543 /*--------------------------------------------------------------------------*/
545 PUBLIC short *make_pair_table(const char *structure)
547 /* returns array representation of structure.
548 table[i] is 0 if unpaired or j if (i.j) pair. */
554 length = (short) strlen(structure);
555 stack = (short *) space(sizeof(short)*(length+1));
556 table = (short *) space(sizeof(short)*(length+2));
559 for (hx=0, i=1; i<=length; i++) {
560 switch (structure[i-1]) {
567 fprintf(stderr, "%s\n", structure);
568 nrerror("unbalanced brackets in make_pair_table");
573 default: /* unpaired base, usually '.' */
579 fprintf(stderr, "%s\n", structure);
580 nrerror("unbalanced brackets in make_pair_table");
586 PUBLIC short *make_pair_table_pk(const char *structure){
593 length = (short) strlen(structure);
594 stack = (short *) space(sizeof(short)*(length+1));
595 stack2 = (short *) space(sizeof(short)*(length+1));
596 table = (short *) space(sizeof(short)*(length+2));
599 for (hx=0, hx2=0, i=1; i<=length; i++) {
600 switch (structure[i-1]) {
607 fprintf(stderr, "%s\n", structure);
608 nrerror("unbalanced '()' brackets in make_pair_table_pk");
619 fprintf(stderr, "%s\n", structure);
620 nrerror("unbalanced '[]' brackets in make_pair_table_pk");
625 default: /* unpaired base, usually '.' */
631 fprintf(stderr, "%s\n", structure);
632 nrerror("unbalanced '()' brackets in make_pair_table_pk");
634 fprintf(stderr, "%s\n", structure);
635 nrerror("unbalanced '[]' brackets in make_pair_table_pk");
642 PUBLIC short *make_pair_table_snoop(const char *structure)
644 /* returns array representation of structure.
645 table[i] is 0 if unpaired or j if (i.j) pair. */
651 length = (short) strlen(structure);
652 stack = (short *) space(sizeof(short)*(length+1));
653 table = (short *) space(sizeof(short)*(length+2));
656 for (hx=0, i=1; i<=length; i++) {
657 switch (structure[i-1]) {
664 fprintf(stderr, "%s\n", structure);
665 nrerror("unbalanced brackets in make_pair_table");
670 default: /* unpaired base, usually '.' */
676 fprintf(stderr, "%s\n", structure);
677 nrerror("unbalanced brackets in make_pair_table");
684 PUBLIC short *alimake_pair_table(const char *structure)
686 /* returns array representation of structure.
687 table[i] is 0 if unpaired or j if (i.j) pair. */
693 length = (short) strlen(structure);
694 stack = (short *) space(sizeof(short)*(length+1));
695 table = (short *) space(sizeof(short)*(length+2));
698 for (hx=0, i=1; i<=length; i++) {
699 switch (structure[i-1]) {
706 fprintf(stderr, "%s\n", structure);
707 nrerror("unbalanced brackets in make_pair_table");
712 default: /* unpaired base, usually '.' */
717 for (hx=0, i=1; i<=length; i++) {
718 switch (structure[i-1]) {
725 fprintf(stderr, "%s\n", structure);
726 nrerror("unbalanced brackets in make_pair_table");
731 default: /* unpaired base, usually '.' */
736 for (hx=0, i=1; i<=length; i++) {
737 switch (structure[i-1]) {
744 fprintf(stderr, "%s\n", structure);
745 nrerror("unbalanced brackets in make_pair_table");
750 default: /* unpaired base, usually '.' */
755 fprintf(stderr, "%s\n", structure);
756 nrerror("unbalanced brackets in make_pair_table");
762 PUBLIC short *copy_pair_table(const short *pt){
763 short *table = (short *)space(sizeof(short) * (pt[0]+2));
764 memcpy(table, pt, sizeof(short)*(pt[0]+2));
769 PUBLIC int *make_loop_index_pt(short *pt){
771 /* number each position by which loop it belongs to (positions start
779 stack = (int *) space(sizeof(int)*(length+1));
780 loop = (int *) space(sizeof(int)*(length+2));
783 for (i=1; i<=length; i++) {
784 if ((pt[i] != 0) && (i < pt[i])) { /* ( */
790 if ((pt[i] != 0) && (i > pt[i])) { /* ) */
793 l = loop[stack[hx-1]]; /* index of enclosing loop */
794 else l=0; /* external loop has index 0 */
796 nrerror("unbalanced brackets in make_pair_table");
805 /*---------------------------------------------------------------------------*/
807 PUBLIC int bp_distance(const char *str1, const char *str2)
809 /* dist = {number of base pairs in one structure but not in the other} */
810 /* same as edit distance with pair_open pair_close as move set */
816 t1 = make_pair_table(str1);
817 t2 = make_pair_table(str2);
819 l = (t1[0]<t2[0])?t1[0]:t2[0]; /* minimum of the two lengths */
831 char *strdup(const char *s) {
834 dup = space(strlen(s)+1);
840 PUBLIC void print_tty_input_seq(void){
841 print_tty_input_seq_str("Input string (upper or lower case)");
844 PUBLIC void print_tty_input_seq_str(const char *s){
845 printf("\n%s; @ to quit\n", s);
846 printf("%s%s\n", scale1, scale2);
847 (void) fflush(stdout);
850 PUBLIC void print_tty_constraint_full(void){
851 print_tty_constraint(VRNA_CONSTRAINT_PIPE | VRNA_CONSTRAINT_DOT | VRNA_CONSTRAINT_X | VRNA_CONSTRAINT_ANG_BRACK | VRNA_CONSTRAINT_RND_BRACK);
854 PUBLIC void print_tty_constraint(unsigned int option){
855 if(!(option & VRNA_CONSTRAINT_NO_HEADER)) printf("Input structure constraints using the following notation:\n");
856 if(option & VRNA_CONSTRAINT_PIPE) printf("| : paired with another base\n");
857 if(option & VRNA_CONSTRAINT_DOT) printf(". : no constraint at all\n");
858 if(option & VRNA_CONSTRAINT_X) printf("x : base must not pair\n");
859 if(option & VRNA_CONSTRAINT_ANG_BRACK) printf("< : base i is paired with a base j<i\n> : base i is paired with a base j>i\n");
860 if(option & VRNA_CONSTRAINT_RND_BRACK) printf("matching brackets ( ): base i pairs base j\n");
863 PUBLIC void str_DNA2RNA(char *sequence){
865 if(sequence != NULL){
866 l = strlen(sequence);
867 for(i = 0; i < l; i++){
868 if(sequence[i] == 'T') sequence[i] = 'U';
869 if(sequence[i] == 't') sequence[i] = 'u';
874 PUBLIC void str_uppercase(char *sequence){
877 l = strlen(sequence);
879 sequence[i] = toupper(sequence[i]);
883 PUBLIC int *get_iindx(unsigned int length){
885 int *idx = (int *)space(sizeof(int) * (length+1));
886 for (i=1; i <= length; i++)
887 idx[i] = (((length + 1 - i) * (length - i))>>1) + length + 1;
891 PUBLIC int *get_indx(unsigned int length){
893 int *idx = (int *)space(sizeof(int) * (length+1));
894 for (i = 1; i <= length; i++)
895 idx[i] = (i*(i-1)) >> 1; /* i(i-1)/2 */
899 PUBLIC void getConstraint(char **cstruc, const char **lines, unsigned int option){
900 int r, i, l, cl, stop;
903 if(option & VRNA_CONSTRAINT_ALL)
904 option |= VRNA_CONSTRAINT_PIPE | VRNA_CONSTRAINT_ANG_BRACK | VRNA_CONSTRAINT_RND_BRACK | VRNA_CONSTRAINT_X | VRNA_CONSTRAINT_G;
906 for(r=i=stop=0;lines[i];i++){
907 l = (int)strlen(lines[i]);
908 c = (char *) space(sizeof(char) * (l+1));
909 (void) sscanf(lines[i], "%s", c);
911 /* line commented out ? */
912 if((*c == '#') || (*c == '%') || (*c == ';') || (*c == '/') || (*c == '*' || (*c == '\0'))){
913 /* skip leading comments only, i.e. do not allow comments inside the constraint */
918 /* check current line for actual constraining structure */
921 case '|': if(!(option & VRNA_CONSTRAINT_PIPE)){
922 warn_user("constraints of type '|' not allowed");
927 case '>': if(!(option & VRNA_CONSTRAINT_ANG_BRACK)){
928 warn_user("constraints of type '<' or '>' not allowed");
933 case ')': if(!(option & VRNA_CONSTRAINT_RND_BRACK)){
934 warn_user("constraints of type '(' or ')' not allowed");
938 case 'x': if(!(option & VRNA_CONSTRAINT_X)){
939 warn_user("constraints of type 'x' not allowed");
943 case '+': if(!(option & VRNA_CONSTRAINT_G)){
944 warn_user("character '+' ignored in structure");
948 case '&': break; /* ignore concatenation char */
949 default: warn_user("unrecognized character in constraint structure");
955 *cstruc = (char *)xrealloc(*cstruc, r*sizeof(char));
956 strcat(*cstruc, ptr);
958 /* stop if not in fasta mode or multiple words on line */
959 if(!(option & VRNA_CONSTRAINT_MULTILINE) || (cl != l)) break;
964 PUBLIC char *extract_record_rest_structure( const char **lines,
966 unsigned int option){
968 char *structure = NULL;
969 int r, i, l, cl, stop;
973 for(r = i = stop = 0; lines[i]; i++){
974 l = (int)strlen(lines[i]);
975 c = (char *) space(sizeof(char) * (l+1));
976 (void) sscanf(lines[i], "%s", c);
979 /* line commented out ? */
980 if((*c == '#') || (*c == '%') || (*c == ';') || (*c == '/') || (*c == '*' || (*c == '\0'))){
981 /* skip leading comments only, i.e. do not allow comments inside the constraint */
986 /* append the structure part to the output */
988 structure = (char *)xrealloc(structure, r*sizeof(char));
989 strcat(structure, c);
991 /* stop if the assumed structure length has been reached */
992 if((length > 0) && (r-1 == length)) break;
993 /* stop if not allowed to read from multiple lines */
994 if(!(option & VRNA_OPTION_MULTILINE)) break;
1002 PUBLIC void constrain_ptypes(const char *constraint, unsigned int length, char *ptype, int *BP, int min_loop_size, unsigned int idx_type){
1008 if(constraint == NULL) return;
1010 n = (int)strlen(constraint);
1012 stack = (int *) space(sizeof(int)*(n+1));
1014 if(!idx_type){ /* index allows access in energy matrices at pos (i,j) via index[j]+i */
1015 index = get_indx(length);
1017 for(hx=0, j=1; j<=n; j++){
1018 switch(constraint[j-1]){
1019 case '|': if(BP) BP[j] = -1;
1021 case 'x': /* can't pair */
1022 for (l=1; l<j-min_loop_size; l++) ptype[index[j]+l] = 0;
1023 for (l=j+min_loop_size+1; l<=(int)length; l++) ptype[index[l]+j] = 0;
1025 case '(': stack[hx++]=j;
1027 case '<': /* pairs upstream */
1028 for (l=1; l<j-min_loop_size; l++) ptype[index[j]+l] = 0;
1030 case ')': if (hx<=0) {
1031 fprintf(stderr, "%s\n", constraint);
1032 nrerror("unbalanced brackets in constraint");
1035 type = ptype[index[j]+i];
1036 for (k=i+1; k<=(int)length; k++) ptype[index[k]+i] = 0;
1037 /* don't allow pairs i<k<j<l */
1038 for (l=j; l<=(int)length; l++)
1039 for (k=i+1; k<=j; k++) ptype[index[l]+k] = 0;
1040 /* don't allow pairs k<i<l<j */
1041 for (l=i; l<=j; l++)
1042 for (k=1; k<=i; k++) ptype[index[l]+k] = 0;
1043 for (k=1; k<j; k++) ptype[index[j]+k] = 0;
1044 ptype[index[j]+i] = (type==0) ? 7 : type;
1046 case '>': /* pairs downstream */
1047 for (l=j+min_loop_size+1; l<=(int)length; l++) ptype[index[l]+j] = 0;
1052 else{ /* index allows access in energy matrices at pos (i,j) via index[i]-j */
1053 index = get_iindx(length);
1055 for(hx=0, j=1; j<=n; j++) {
1056 switch (constraint[j-1]) {
1057 case 'x': /* can't pair */
1058 for (l=1; l<j-min_loop_size; l++) ptype[index[l]-j] = 0;
1059 for (l=j+min_loop_size+1; l<=(int)length; l++) ptype[index[j]-l] = 0;
1061 case '(': stack[hx++]=j;
1063 case '<': /* pairs upstream */
1064 for (l=1; l<j-min_loop_size; l++) ptype[index[l]-j] = 0;
1066 case ')': if (hx<=0) {
1067 fprintf(stderr, "%s\n", constraint);
1068 nrerror("unbalanced brackets in constraints");
1071 type = ptype[index[i]-j];
1072 /* don't allow pairs i<k<j<l */
1073 for (k=i; k<=j; k++)
1074 for (l=j; l<=(int)length; l++) ptype[index[k]-l] = 0;
1075 /* don't allow pairs k<i<l<j */
1076 for (k=1; k<=i; k++)
1077 for (l=i; l<=j; l++) ptype[index[k]-l] = 0;
1078 ptype[index[i]-j] = (type==0) ? 7 : type;
1080 case '>': /* pairs downstream */
1081 for (l=j+min_loop_size+1; l<=(int)length; l++) ptype[index[j]-l] = 0;
1087 fprintf(stderr, "%s\n", constraint);
1088 nrerror("unbalanced brackets in constraint string");
1094 /* get a matrix containing the number of basepairs of a reference structure for each interval [i,j] with i<j
1095 * access it via iindx!!!
1097 PUBLIC unsigned int *make_referenceBP_array(short *reference_pt, unsigned int turn){
1098 unsigned int i,j,k,ij,length;
1100 unsigned int *array;
1102 length = (unsigned int)reference_pt[0];
1103 size = ((length+1)*(length+2))/2;
1104 iindx = get_iindx(length);
1105 array = (unsigned int *) space(sizeof(unsigned int)*size); /* matrix containing number of basepairs of reference structure1 in interval [i,j] */;
1106 for (k=0; k<=turn; k++)
1107 for (i=1; i<=length-k; i++) {
1113 for (i = length-turn-1; i >= 1; i--)
1114 for (j = i+turn+1; j <= length; j++){
1118 if((i<=(unsigned int)reference_pt[j]) && ((unsigned int)reference_pt[j] < j))
1126 PUBLIC unsigned int *compute_BPdifferences(short *pt1, short *pt2, unsigned int turn){
1127 unsigned int *array;
1128 unsigned int n, size, i, j, ij, d;
1129 n = (unsigned int)pt1[0];
1130 size = ((n+1)*(n+2))/2;
1131 array = (unsigned int *)space(sizeof(unsigned int) * size);
1132 int *iindx = get_iindx(n);
1133 for(i = n - turn - 1; i>=1; i--){
1135 for(j = i+turn+1; j <= n; j++){
1138 if(pt1[j] != pt2[j]){
1139 if(i <= (unsigned int)pt1[j] && (unsigned int)pt1[j] < j){
1140 /* we got an additional base pair in reference structure 1 */
1143 if(i <= (unsigned int)pt2[j] && (unsigned int)pt2[j] < j){
1144 /* we got another base pair in reference structure 2 */