binaries/src/ViennaRNA/lib/utils.c

   1 /*
   2                                utils.c
   3
   4                  c  Ivo L Hofacker and Walter Fontana
   5                           Vienna RNA package
   6 */
   7 /* Last changed Time-stamp: <2008-11-25 16:34:36 ivo> */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <ctype.h>
  12 #include <errno.h>
  13 #include <time.h>
  14 #include <string.h>
  15 #include "../config.h"
  16 #include "utils.h"
  17
  18 #ifdef WITH_DMALLOC
  19 #include "dmalloc.h"
  20 #endif
  21
  22 #define PRIVATE  static
  23 #define PUBLIC
  24
  25 /*@notnull@ @only@*/
  26 PUBLIC unsigned short xsubi[3];
  27
  28 PRIVATE char  scale1[] = "....,....1....,....2....,....3....,....4";
  29 PRIVATE char  scale2[] = "....,....5....,....6....,....7....,....8";
  30 PRIVATE char  *inbuf = NULL;
  31
  32 PRIVATE char  *inbuf2 = NULL;
  33 PRIVATE unsigned int typebuf2 = 0;
  34
  35 /*-------------------------------------------------------------------------*/
  36
  37 PUBLIC void *space(unsigned size) {
  38   void *pointer;
  39
  40   if ( (pointer = (void *) calloc(1, (size_t) size)) == NULL) {
  41 #ifdef EINVAL
  42     if (errno==EINVAL) {
  43       fprintf(stderr,"SPACE: requested size: %d\n", size);
  44       nrerror("SPACE allocation failure -> EINVAL");
  45     }
  46     if (errno==ENOMEM)
  47 #endif
  48       nrerror("SPACE allocation failure -> no memory");
  49   }
  50   return  pointer;
  51 }
  52
  53 #ifdef WITH_DMALLOC
  54 #define space(S) calloc(1,(S))
  55 #endif
  56
  57 #undef xrealloc
  58 /* dmalloc.h #define's xrealloc */
  59 void *xrealloc (void *p, unsigned size) {
  60   if (p == 0)
  61     return space(size);
  62   p = (void *) realloc(p, size);
  63   if (p == NULL) {
  64 #ifdef EINVAL
  65     if (errno==EINVAL) {
  66       fprintf(stderr,"xrealloc: requested size: %d\n", size);
  67       nrerror("xrealloc allocation failure -> EINVAL");
  68     }
  69     if (errno==ENOMEM)
  70 #endif
  71       nrerror("xrealloc allocation failure -> no memory");
  72   }
  73   return p;
  74 }
  75
  76 /*------------------------------------------------------------------------*/
  77
  78 PUBLIC void nrerror(const char message[])       /* output message upon error */
  79 {
  80   fprintf(stderr, "ERROR: %s\n", message);
  81   exit(EXIT_FAILURE);
  82 }
  83
  84 PUBLIC void warn_user(const char message[]){
  85   fprintf(stderr, "WARNING: %s\n", message);
  86 }
  87
  88 /*------------------------------------------------------------------------*/
  89 PUBLIC void init_rand(void)
  90 {
  91   time_t t;
  92   (void) time(&t);
  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);
  96 #ifndef HAVE_ERAND48
  97   srand((unsigned int) t);
  98 #endif
  99 }
 100
 101 /*------------------------------------------------------------------------*/
 102
 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 */
 107 {
 108 #ifdef HAVE_ERAND48
 109   extern double erand48(unsigned short[]);
 110   return erand48(xsubi);
 111 #else
 112   return ((double) rand())/RAND_MAX;
 113 #endif
 114 }
 115
 116 /*------------------------------------------------------------------------*/
 117
 118 PUBLIC int int_urn(int from, int to)
 119 {
 120   return ( ( (int) (urn()*(to-from+1)) ) + from );
 121 }
 122
 123 /*------------------------------------------------------------------------*/
 124
 125 PUBLIC void filecopy(FILE *from, FILE *to)
 126 {
 127   int c;
 128
 129   while ((c = getc(from)) != EOF) (void)putc(c, to);
 130 }
 131
 132 /*-----------------------------------------------------------------*/
 133
 134 PUBLIC char *time_stamp(void)
 135 {
 136   time_t  cal_time;
 137
 138   cal_time = time(NULL);
 139   return ( ctime(&cal_time) );
 140 }
 141
 142 /*-----------------------------------------------------------------*/
 143
 144 PUBLIC char *random_string(int l, const char symbols[])
 145 {
 146   char *r;
 147   int   i, rn, base;
 148
 149   base = (int) strlen(symbols);
 150   r = (char *) space(sizeof(char)*(l+1));
 151
 152   for (i = 0; i < l; i++) {
 153     rn = (int) (urn()*base);  /* [0, base-1] */
 154     r[i] = symbols[rn];
 155   }
 156   r[l] = '\0';
 157   return r;
 158 }
 159
 160 /*-----------------------------------------------------------------*/
 161
 162 PUBLIC int   hamming(const char *s1, const char *s2)
 163 {
 164   int h=0;
 165
 166   for (; *s1 && *s2; s1++, s2++)
 167     if (*s1 != *s2) h++;
 168   return h;
 169 }
 170
 171 PUBLIC int   hamming_bound(const char *s1, const char *s2, int boundary)
 172 {
 173   int h=0;
 174
 175   for (; *s1 && *s2 && boundary; s1++, s2++, boundary--)
 176     if (*s1 != *s2) h++;
 177   return h;
 178 }
 179 /*-----------------------------------------------------------------*/
 180
 181 PUBLIC char *get_line(FILE *fp) /* reads lines of arbitrary length from fp */
 182 {
 183   char s[512], *line, *cp;
 184   int len=0, size=0, l;
 185   line=NULL;
 186   do {
 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);
 191     if (l+1>size) {
 192       size = (int)((l+1)*1.2);
 193       line = (char *) xrealloc(line, size*sizeof(char));
 194     }
 195     strcat(line+len, s);
 196     len=l;
 197   } while(cp==NULL);
 198
 199   return line;
 200 }
 201
 202 PUBLIC int  skip_comment_lines(char **line){
 203   if((*line = get_line(stdin))==NULL) return -1;
 204
 205   while((**line=='*')||(**line=='\0')){
 206     free(*line);
 207     if((*line = get_line(stdin))==NULL) return -1;
 208   }
 209   return 0;
 210 }
 211
 212 PUBLIC  unsigned int get_input_line(char **string, unsigned int option){
 213   char  *line;
 214   int   i, l, r;
 215
 216   /*
 217   * read lines until informative data appears or
 218   * report an error if anything goes wrong
 219   */
 220   if((line = get_line(stdin))==NULL) return VRNA_INPUT_ERROR;
 221
 222   if(!(option & VRNA_INPUT_NOSKIP_COMMENTS))
 223     while((*line=='*')||(*line=='\0')){
 224       free(line);
 225       if((line = get_line(stdin))==NULL) return VRNA_INPUT_ERROR;
 226     }
 227
 228   l = (int) strlen(line);
 229
 230   /* break on '@' sign if not disabled */
 231   if(*line == '@'){
 232     free(line);
 233     return VRNA_INPUT_QUIT;
 234   }
 235   /* print line read if not disabled */
 236   /* if(!(option & VRNA_INPUT_NOPRINT)) printf("%s\n", line); */
 237
 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;
 243       else                      break;
 244     }
 245     line[(i >= 0) ? (i+1) : 0] = '\0';
 246   }
 247
 248   if(*line == '>'){
 249     /* fasta header */
 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);
 254     if(i > 0){
 255       i       = (int)     strlen(*string);
 256       *string = (char *)  xrealloc(*string, (i+1)*sizeof(char));
 257       free(line);
 258       return r;
 259     }
 260     else{
 261       free(line);
 262       free(*string);
 263       *string = NULL;
 264       return VRNA_INPUT_ERROR;
 265     }
 266   }
 267   else{
 268     *string = strdup(line);
 269     free(line);
 270   }
 271   return VRNA_INPUT_MISC;
 272 }
 273
 274 PUBLIC  unsigned int get_multi_input_line(char **string, unsigned int option){
 275   char  *line;
 276   int   i, l;
 277   int   state = 0;
 278   int   str_length = 0;
 279
 280   line = (inbuf) ? inbuf : get_line(stdin);
 281   inbuf = NULL;
 282   do{
 283
 284     /*
 285     * read lines until informative data appears or
 286     * report an error if anything goes wrong
 287     */
 288     if(!line) return VRNA_INPUT_ERROR;
 289
 290     l = (int)strlen(line);
 291
 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;
 297         else                      break;
 298       }
 299       line[(i >= 0) ? (i+1) : 0] = '\0';
 300     }
 301
 302     l           = (int)strlen(line);
 303     str_length  = (*string) ? (int) strlen(*string) : 0;
 304
 305     switch(*line){
 306       case  '@':    /* user abort */
 307                     if(state) inbuf = line;
 308                     else      free(line);
 309                     return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_QUIT;
 310
 311       case  '\0':   /* empty line */
 312                     if(option & VRNA_INPUT_NOSKIP_BLANK_LINES){
 313                       if(state) inbuf = line;
 314                       else      free(line);
 315                       return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_BLANK_LINE;
 316                     }
 317                     break;
 318
 319       case  '#': case  '%': case  ';': case  '/': case  '*': case ' ':
 320                     /* comments */
 321                     if(option & VRNA_INPUT_NOSKIP_COMMENTS){
 322                       if(state) inbuf   = line;
 323                       else      *string = line;
 324                       return (state == 2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_COMMENT;
 325                     }
 326                     break;
 327
 328       case  '>':    /* fasta header */
 329                     if(state) inbuf   = line;
 330                     else      *string = line;
 331                     return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_FASTA_HEADER;
 332
 333       case  'x': case 'e': case 'l': case '&':   /* seems to be a constraint or line starting with second sequence for dimer calculations */
 334                     i = 1;
 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 */
 338
 339                     if(
 340                             ((line[i]>64) && (line[i]<91))  /* A-Z */
 341                         ||  ((line[i]>96) && (line[i]<123)) /* a-z */
 342                       ){
 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;}
 346                         else{
 347                           *string = (char *)xrealloc(*string, sizeof(char) * (str_length + l + 1));
 348                           strcpy(*string + str_length, line);
 349                           state = 1;
 350                         }
 351                         break;
 352                       }
 353                       /* otherwise return line read */
 354                       else{ *string = line; return VRNA_INPUT_SEQUENCE;}
 355                     }
 356                     /* mmmh? it really seems to be a constraint */
 357                     /* fallthrough */
 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){
 362                       if(state == 1){
 363                         inbuf = line;
 364                         return VRNA_INPUT_SEQUENCE;
 365                       }
 366                       else{
 367                         *string = (char *)xrealloc(*string, sizeof(char) * (str_length + l + 1));
 368                         strcpy(*string + str_length, line);
 369                         state = 2;
 370                       }
 371                     }
 372                     /* or we return it as it is */
 373                     else{
 374                       *string = line;
 375                       return VRNA_INPUT_CONSTRAINT;
 376                     }
 377                     break;
 378       default:      if(option & VRNA_INPUT_FASTA_HEADER){
 379                       /* are we already in sequence mode? */
 380                       if(state == 2){
 381                         inbuf = line;
 382                         return VRNA_INPUT_CONSTRAINT;
 383                       }
 384                       else{
 385                         *string = (char *)xrealloc(*string, sizeof(char) * (str_length + l + 1));
 386                         strcpy(*string + str_length, line);
 387                         state = 1;
 388                       }
 389                     }
 390                     /* otherwise return line read */
 391                     else{
 392                       *string = line;
 393                       return VRNA_INPUT_SEQUENCE;
 394                     }
 395     }
 396     free(line);
 397     line = get_line(stdin);
 398   }while(line);
 399
 400   return (state==2) ? VRNA_INPUT_CONSTRAINT : (state==1) ? VRNA_INPUT_SEQUENCE : VRNA_INPUT_ERROR;
 401 }
 402
 403 PUBLIC  unsigned int read_record(char **header, char **sequence, char ***rest, unsigned int options){
 404   unsigned int  input_type, return_type, tmp_type;
 405   int           rest_count;
 406   char          *input_string;
 407
 408   rest_count    = 0;
 409   return_type   = tmp_type = 0;
 410   input_string  = *header = *sequence = NULL;
 411   *rest         = (char **)space(sizeof(char *));
 412
 413   /* remove unnecessary option flags from options variable... */
 414   options &= ~VRNA_INPUT_FASTA_HEADER;
 415
 416   /* read first input or last buffered input */
 417   if(typebuf2){
 418     input_type    = typebuf2;
 419     input_string  = inbuf2;
 420     typebuf2      = 0;
 421     inbuf2        = NULL;
 422   }
 423   else input_type  = get_multi_input_line(&input_string, options);
 424
 425   if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) return input_type;
 426
 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;
 432   }
 433
 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;
 437     input_string  = NULL;
 438     /* get next data-block with fasta support if not explicitely forbidden by VRNA_INPUT_NO_SPAN */
 439     input_type  = get_multi_input_line(
 440                     &input_string,
 441                     ((options & VRNA_INPUT_NO_SPAN) ? 0 : VRNA_INPUT_FASTA_HEADER) | options
 442                   );
 443     if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) return (return_type | input_type);
 444   }
 445
 446   if(input_type & VRNA_INPUT_SEQUENCE){
 447     return_type  |= VRNA_INPUT_SEQUENCE; /* remember that we've read a sequence */
 448     *sequence     = input_string;
 449     input_string  = NULL;
 450   } else nrerror("sequence input missing");
 451
 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;
 460       input_string = NULL;
 461     }
 462     /*
 463     if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) return input_type;
 464     */
 465
 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
 469     */
 470     inbuf2 = input_string;
 471     typebuf2 = input_type;
 472   }
 473   (*rest)[rest_count] = NULL;
 474   return (return_type);
 475 }
 476
 477
 478 /*-----------------------------------------------------------------*/
 479
 480 PUBLIC char *pack_structure(const char *struc) {
 481   /* 5:1 compression using base 3 encoding */
 482   int i,j,l,pi;
 483   unsigned char *packed;
 484
 485   l = (int) strlen(struc);
 486   packed = (unsigned char *) space(((l+4)/5+1)*sizeof(unsigned char));
 487
 488   j=i=pi=0;
 489   while (i<l) {
 490     register int p;
 491     for (p=pi=0; pi<5; pi++) {
 492       p *= 3;
 493       switch (struc[i]) {
 494       case '(':
 495       case '\0':
 496         break;
 497       case '.':
 498         p++;
 499         break;
 500       case ')':
 501         p += 2;
 502         break;
 503       default: nrerror("pack_structure: illegal charcter in structure");
 504       }
 505       if (i<l) i++;
 506     }
 507     packed[j++] = (unsigned char) (p+1); /* never use 0, so we can use
 508                                             strcmp()  etc. */
 509   }
 510   packed[j] = '\0';      /* for str*() functions */
 511   return (char *) packed;
 512 }
 513
 514 PUBLIC char *unpack_structure(const char *packed) {
 515   /* 5:1 compression using base 3 encoding */
 516   int i,j,l;
 517   char *struc;
 518   unsigned const char *pp;
 519   char code[3] = {'(', '.', ')'};
 520
 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 */
 524
 525   for (i=j=0; i<l; i++) {
 526     register int p, c, k;
 527
 528     p = (int) pp[i] - 1;
 529     for (k=4; k>=0; k--) {
 530       c = p % 3;
 531       p /= 3;
 532       struc[j+k] = code[c];
 533     }
 534     j += 5;
 535   }
 536   struc[j--] = '\0';
 537   while (struc[j] == '(') /* strip trailing ( */
 538     struc[j--] = '\0';
 539
 540   return struc;
 541 }
 542
 543 /*--------------------------------------------------------------------------*/
 544
 545 PUBLIC short *make_pair_table(const char *structure)
 546 {
 547     /* returns array representation of structure.
 548        table[i] is 0 if unpaired or j if (i.j) pair.  */
 549    short i,j,hx;
 550    short length;
 551    short *stack;
 552    short *table;
 553
 554    length = (short) strlen(structure);
 555    stack = (short *) space(sizeof(short)*(length+1));
 556    table = (short *) space(sizeof(short)*(length+2));
 557    table[0] = length;
 558
 559    for (hx=0, i=1; i<=length; i++) {
 560       switch (structure[i-1]) {
 561        case '(':
 562          stack[hx++]=i;
 563          break;
 564        case ')':
 565          j = stack[--hx];
 566          if (hx<0) {
 567             fprintf(stderr, "%s\n", structure);
 568             nrerror("unbalanced brackets in make_pair_table");
 569          }
 570          table[i]=j;
 571          table[j]=i;
 572          break;
 573        default:   /* unpaired base, usually '.' */
 574          table[i]= 0;
 575          break;
 576       }
 577    }
 578    if (hx!=0) {
 579       fprintf(stderr, "%s\n", structure);
 580       nrerror("unbalanced brackets in make_pair_table");
 581    }
 582    free(stack);
 583    return(table);
 584 }
 585
 586 PUBLIC short *make_pair_table_pk(const char *structure){
 587    short i,j,hx, hx2;
 588    short length;
 589    short *stack;
 590    short *stack2;
 591    short *table;
 592
 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));
 597    table[0] = length;
 598
 599    for (hx=0, hx2=0, i=1; i<=length; i++) {
 600       switch (structure[i-1]) {
 601        case '(':
 602          stack[hx++]=i;
 603          break;
 604        case ')':
 605          j = stack[--hx];
 606          if (hx<0) {
 607             fprintf(stderr, "%s\n", structure);
 608             nrerror("unbalanced '()' brackets in make_pair_table_pk");
 609          }
 610          table[i]=j;
 611          table[j]=i;
 612          break;
 613        case '[':
 614          stack2[hx2++]=i;
 615          break;
 616        case ']':
 617          j = stack2[--hx2];
 618          if (hx2<0) {
 619             fprintf(stderr, "%s\n", structure);
 620             nrerror("unbalanced '[]' brackets in make_pair_table_pk");
 621          }
 622          table[i]=j;
 623          table[j]=i;
 624          break;
 625        default:   /* unpaired base, usually '.' */
 626          table[i]= 0;
 627          break;
 628       }
 629    }
 630    if (hx!=0) {
 631       fprintf(stderr, "%s\n", structure);
 632       nrerror("unbalanced '()' brackets in make_pair_table_pk");
 633    } else if (hx2!=0) {
 634       fprintf(stderr, "%s\n", structure);
 635       nrerror("unbalanced '[]' brackets in make_pair_table_pk");
 636    }
 637    free(stack);
 638    free(stack2);
 639    return(table);
 640 }
 641
 642 PUBLIC short *make_pair_table_snoop(const char *structure)
 643 {
 644     /* returns array representation of structure.
 645        table[i] is 0 if unpaired or j if (i.j) pair.  */
 646    short i,j,hx;
 647    short length;
 648    short *stack;
 649    short *table;
 650
 651    length = (short) strlen(structure);
 652    stack = (short *) space(sizeof(short)*(length+1));
 653    table = (short *) space(sizeof(short)*(length+2));
 654    table[0] = length;
 655
 656    for (hx=0, i=1; i<=length; i++) {
 657      switch (structure[i-1]) {
 658      case '<':
 659        stack[hx++]=i;
 660        break;
 661      case '>':
 662        j = stack[--hx];
 663        if (hx<0) {
 664          fprintf(stderr, "%s\n", structure);
 665          nrerror("unbalanced brackets in make_pair_table");
 666        }
 667        table[i]=j;
 668        table[j]=i;
 669        break;
 670      default:   /* unpaired base, usually '.' */
 671        table[i]= table[i];
 672        break;
 673      }
 674    }
 675    if (hx!=0) {
 676      fprintf(stderr, "%s\n", structure);
 677      nrerror("unbalanced brackets in make_pair_table");
 678    }
 679    free(stack);
 680    return table ;
 681 }
 682
 683
 684 PUBLIC short *alimake_pair_table(const char *structure)
 685 {
 686     /* returns array representation of structure.
 687        table[i] is 0 if unpaired or j if (i.j) pair.  */
 688    short i,j,hx;
 689    short length;
 690    short *stack;
 691    short *table;
 692
 693    length = (short) strlen(structure);
 694    stack = (short *) space(sizeof(short)*(length+1));
 695    table = (short *) space(sizeof(short)*(length+2));
 696    table[0] = length;
 697
 698    for (hx=0, i=1; i<=length; i++) {
 699       switch (structure[i-1]) {
 700        case '(':
 701          stack[hx++]=i;
 702          break;
 703        case ')':
 704          j = stack[--hx];
 705          if (hx<0) {
 706             fprintf(stderr, "%s\n", structure);
 707             nrerror("unbalanced brackets in make_pair_table");
 708          }
 709          table[i]=j;
 710          table[j]=i;
 711          break;
 712        default:   /* unpaired base, usually '.' */
 713          table[i]= 0;
 714          break;
 715       }
 716    }
 717    for (hx=0, i=1; i<=length; i++) {
 718       switch (structure[i-1]) {
 719        case '<':
 720          stack[hx++]=i;
 721          break;
 722        case '>':
 723          j = stack[--hx];
 724          if (hx<0) {
 725             fprintf(stderr, "%s\n", structure);
 726             nrerror("unbalanced brackets in make_pair_table");
 727          }
 728          table[i]=j;
 729          table[j]=i;
 730          break;
 731        default:   /* unpaired base, usually '.' */
 732          table[i]= table[i];
 733          break;
 734       }
 735    }
 736    for (hx=0, i=1; i<=length; i++) {
 737      switch (structure[i-1]) {
 738      case '[':
 739        stack[hx++]=i;
 740        break;
 741      case ']':
 742        j = stack[--hx];
 743        if (hx<0) {
 744          fprintf(stderr, "%s\n", structure);
 745          nrerror("unbalanced brackets in make_pair_table");
 746        }
 747        table[i]=j;
 748        table[j]=i;
 749        break;
 750      default:   /* unpaired base, usually '.' */
 751        break;
 752      }
 753    }
 754    if (hx!=0) {
 755       fprintf(stderr, "%s\n", structure);
 756       nrerror("unbalanced brackets in make_pair_table");
 757    }
 758    free(stack);
 759    return(table);
 760 }
 761
 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));
 765   return table;
 766 }
 767
 768
 769 PUBLIC int *make_loop_index_pt(short *pt){
 770
 771   /* number each position by which loop it belongs to (positions start
 772      at 1) */
 773   int i,hx,l,nl;
 774   int length;
 775   int *stack = NULL;
 776   int *loop = NULL;
 777
 778   length = pt[0];
 779   stack  = (int *) space(sizeof(int)*(length+1));
 780   loop   = (int *) space(sizeof(int)*(length+2));
 781   hx=l=nl=0;
 782
 783   for (i=1; i<=length; i++) {
 784     if ((pt[i] != 0) && (i < pt[i])) { /* ( */
 785       nl++; l=nl;
 786       stack[hx++]=i;
 787     }
 788     loop[i]=l;
 789
 790     if ((pt[i] != 0) && (i > pt[i])) { /* ) */
 791       --hx;
 792       if (hx>0)
 793         l = loop[stack[hx-1]];  /* index of enclosing loop   */
 794       else l=0;                 /* external loop has index 0 */
 795       if (hx<0) {
 796         nrerror("unbalanced brackets in make_pair_table");
 797       }
 798     }
 799   }
 800   loop[0] = nl;
 801   free(stack);
 802   return (loop);
 803 }
 804
 805 /*---------------------------------------------------------------------------*/
 806
 807 PUBLIC int bp_distance(const char *str1, const char *str2)
 808 {
 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 */
 811    int dist;
 812    short i,l;
 813    short *t1, *t2;
 814
 815    dist = 0;
 816    t1 = make_pair_table(str1);
 817    t2 = make_pair_table(str2);
 818
 819    l = (t1[0]<t2[0])?t1[0]:t2[0];    /* minimum of the two lengths */
 820
 821    for (i=1; i<=l; i++)
 822      if (t1[i]!=t2[i]) {
 823        if (t1[i]>i) dist++;
 824        if (t2[i]>i) dist++;
 825      }
 826    free(t1); free(t2);
 827    return dist;
 828 }
 829
 830 #ifndef HAVE_STRDUP
 831 char *strdup(const char *s) {
 832   char *dup;
 833
 834   dup = space(strlen(s)+1);
 835   strcpy(dup, s);
 836   return(dup);
 837 }
 838 #endif
 839
 840 PUBLIC  void  print_tty_input_seq(void){
 841   print_tty_input_seq_str("Input string (upper or lower case)");
 842 }
 843
 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);
 848 }
 849
 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);
 852 }
 853
 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");
 861 }
 862
 863 PUBLIC  void  str_DNA2RNA(char *sequence){
 864   unsigned int l, i;
 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';
 870     }
 871   }
 872 }
 873
 874 PUBLIC void str_uppercase(char *sequence){
 875   unsigned int l, i;
 876   if(sequence){
 877     l = strlen(sequence);
 878     for(i=0;i<l;i++)
 879       sequence[i] = toupper(sequence[i]);
 880   }
 881 }
 882
 883 PUBLIC int *get_iindx(unsigned int length){
 884   int i;
 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;
 888   return idx;
 889 }
 890
 891 PUBLIC int *get_indx(unsigned int length){
 892   unsigned int i;
 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 */
 896   return idx;
 897 }
 898
 899 PUBLIC void getConstraint(char **cstruc, const char **lines, unsigned int option){
 900   int r, i, l, cl, stop;
 901   char *c, *ptr;
 902   if(lines){
 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;
 905
 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);
 910       cl  = (int)strlen(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 */
 914         if(!r)  continue;
 915         else    break;
 916       }
 917
 918       /* check current line for actual constraining structure */
 919       for(ptr = c;*c;c++){
 920         switch(*c){
 921           case '|':   if(!(option & VRNA_CONSTRAINT_PIPE)){
 922                         warn_user("constraints of type '|' not allowed");
 923                         *c = '.';
 924                       }
 925                       break;
 926           case '<':
 927           case '>':   if(!(option & VRNA_CONSTRAINT_ANG_BRACK)){
 928                         warn_user("constraints of type '<' or '>' not allowed");
 929                         *c = '.';
 930                       }
 931                       break;
 932           case '(':
 933           case ')':   if(!(option & VRNA_CONSTRAINT_RND_BRACK)){
 934                         warn_user("constraints of type '(' or ')' not allowed");
 935                         *c = '.';
 936                       }
 937                       break;
 938           case 'x':   if(!(option & VRNA_CONSTRAINT_X)){
 939                         warn_user("constraints of type 'x' not allowed");
 940                         *c = '.';
 941                       }
 942                       break;
 943           case '+':   if(!(option & VRNA_CONSTRAINT_G)){
 944                         warn_user("character '+' ignored in structure");
 945                         *c = '.';
 946                       }
 947           case '.':   break;
 948           case '&':   break; /* ignore concatenation char */
 949           default:    warn_user("unrecognized character in constraint structure");
 950                       break;
 951         }
 952       }
 953
 954       r += cl+1;
 955       *cstruc = (char *)xrealloc(*cstruc, r*sizeof(char));
 956       strcat(*cstruc, ptr);
 957       free(ptr);
 958       /* stop if not in fasta mode or multiple words on line */
 959       if(!(option & VRNA_CONSTRAINT_MULTILINE) || (cl != l)) break;
 960     }
 961   }
 962 }
 963
 964 PUBLIC char *extract_record_rest_structure( const char **lines,
 965                                             unsigned int length,
 966                                             unsigned int option){
 967
 968   char *structure = NULL;
 969   int r, i, l, cl, stop;
 970   char *c;
 971
 972   if(lines){
 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);
 977       cl  = (int)strlen(c);
 978
 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 */
 982         if(!r)  continue;
 983         else    break;
 984       }
 985
 986       /* append the structure part to the output */
 987       r += cl+1;
 988       structure = (char *)xrealloc(structure, r*sizeof(char));
 989       strcat(structure, c);
 990       free(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;
 995     }
 996   }
 997   return structure;
 998 }
 999
1000
1001
1002 PUBLIC void constrain_ptypes(const char *constraint, unsigned int length, char *ptype, int *BP, int min_loop_size, unsigned int idx_type){
1003   int n,i,j,k,l;
1004   int hx, *stack;
1005   char type;
1006   int *index;
1007
1008   if(constraint == NULL) return;
1009
1010   n = (int)strlen(constraint);
1011
1012   stack = (int *) space(sizeof(int)*(n+1));
1013
1014   if(!idx_type){ /* index allows access in energy matrices at pos (i,j) via index[j]+i */
1015     index = get_indx(length);
1016
1017     for(hx=0, j=1; j<=n; j++){
1018       switch(constraint[j-1]){
1019         case '|':   if(BP) BP[j] = -1;
1020                     break;
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;
1024                     break;
1025         case '(':   stack[hx++]=j;
1026                     /* fallthrough */
1027         case '<':   /* pairs upstream */
1028                     for (l=1; l<j-min_loop_size; l++) ptype[index[j]+l] = 0;
1029                     break;
1030         case ')':   if (hx<=0) {
1031                       fprintf(stderr, "%s\n", constraint);
1032                       nrerror("unbalanced brackets in constraint");
1033                     }
1034                     i = stack[--hx];
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;
1045                     /* fallthrough */
1046         case '>':   /* pairs downstream */
1047                     for (l=j+min_loop_size+1; l<=(int)length; l++) ptype[index[l]+j] = 0;
1048                     break;
1049       }
1050     }
1051   }
1052   else{ /* index allows access in energy matrices at pos (i,j) via index[i]-j */
1053     index = get_iindx(length);
1054
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;
1060                     break;
1061         case '(':   stack[hx++]=j;
1062                     /* fallthrough */
1063         case '<':   /* pairs upstream */
1064                     for (l=1; l<j-min_loop_size; l++) ptype[index[l]-j] = 0;
1065                     break;
1066         case ')':   if (hx<=0) {
1067                       fprintf(stderr, "%s\n", constraint);
1068                       nrerror("unbalanced brackets in constraints");
1069                     }
1070                     i = stack[--hx];
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;
1079                     /* fallthrough */
1080         case '>':   /* pairs downstream */
1081                     for (l=j+min_loop_size+1; l<=(int)length; l++) ptype[index[j]-l] = 0;
1082                     break;
1083       }
1084     }
1085   }
1086   if (hx!=0) {
1087     fprintf(stderr, "%s\n", constraint);
1088     nrerror("unbalanced brackets in constraint string");
1089   }
1090   free(index);
1091   free(stack);
1092 }
1093
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!!!
1096 */
1097 PUBLIC unsigned int *make_referenceBP_array(short *reference_pt, unsigned int turn){
1098   unsigned int i,j,k,ij,length;
1099   int *iindx;
1100   unsigned int *array;
1101   unsigned int size;
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++) {
1108       j=i+k;
1109       ij = iindx[i]-j;
1110       array[ij] = 0;
1111     }
1112
1113   for (i = length-turn-1; i >= 1; i--)
1114     for (j = i+turn+1; j <= length; j++){
1115       int bps;
1116       ij = iindx[i]-j;
1117       bps = array[ij+1];
1118       if((i<=(unsigned int)reference_pt[j]) && ((unsigned int)reference_pt[j] < j))
1119         bps++;
1120       array[ij] = bps;
1121     }
1122   free(iindx);
1123   return array;
1124 }
1125
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--){
1134     d = 0;
1135     for(j = i+turn+1; j <= n; j++){
1136       ij = iindx[i]-j;
1137       d = array[ij+1];
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 */
1141           d++;
1142         }
1143         if(i <= (unsigned int)pt2[j] && (unsigned int)pt2[j] < j){
1144           /* we got another base pair in reference structure 2 */
1145           d++;
1146         }
1147       }
1148       array[ij] = d;
1149
1150     }
1151   }
1152   free(iindx);
1153   return array;
1154 }