binaries/src/tcoffee/t_coffee_source/util.c

   1 #define FILE_CHECK 1
   2 #include <stddef.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <ctype.h>
   6 #include <string.h>
   7 #include <limits.h>
   8 #include <float.h>
   9 #include <math.h>
  10 #include <stdarg.h>
  11 #include <fcntl.h>
  12 #include <sys/file.h>
  13
  14
  15 #include "io_lib_header.h"
  16 #include "util_lib_header.h"
  17 #include "define_header.h"
  18 #include "perl_header_lib.h"
  19
  20 int my_vsscanf(char *buf, char *fmt, va_list parms);
  21 static int get_vtmpnam2_root();
  22
  23 static char CACHE_4_TCOFFEE[1000];
  24 static char TMP_4_TCOFFEE[1000];
  25 static char DIR_4_TCOFFEE[1000];
  26 static int global_exit_signal;
  27 static int no_error_report;
  28 /*********************************************************************/
  29 /*                                                                   */
  30 /*                                  DICHOTOMY                        */
  31 /*                                                                   */
  32 /*                                                                   */
  33 /*********************************************************************/
  34 double dichotomy (double value, double target_value, double middle, double *bottom,double *top)
  35 {
  36   if ( value> target_value)top[0]=middle;
  37   else if ( value<target_value)bottom[0]=middle;
  38   return (top[0]-bottom[0])/2+bottom[0];
  39 }
  40
  41
  42 /*********************************************************************/
  43 /*                                                                   */
  44 /*                                   QSORT                           */
  45 /*                                                                   */
  46 /*                                                                   */
  47 /*********************************************************************/
  48
  49 #define MAXSTACK (sizeof(size_t) * CHAR_BIT)
  50 static void exchange(void *a, void *b, size_t size) {
  51     size_t i;
  52     int *ia;
  53     int *ib;
  54     char *ca;
  55     char *cb;
  56
  57     /******************
  58      *  exchange a,b  *
  59      ******************/
  60     ia=(int*)a;
  61     ib=(int*)b;
  62     for (i = sizeof(int); i <= size; i += sizeof(int)) {
  63         int t =ia[0];
  64         ia[0]=ib[0];
  65         ib[0]=t;
  66         ia++; ib++;
  67
  68     }
  69     ca=(char*)ia;
  70     cb=(char*)ib;
  71     for (i = i - sizeof(int) + 1; i <= size; i++) {
  72
  73         char t=ca[0];
  74         ca[0]=cb[0];
  75         cb[0]=t;
  76         ca++;cb++;
  77
  78     }
  79 }
  80
  81 void qsort(void *base, size_t nmemb, size_t size,
  82         int (*compar)(const void *, const void *)) {
  83     void *lbStack[MAXSTACK], *ubStack[MAXSTACK];
  84     int sp;
  85     unsigned int offset;
  86
  87     /********************
  88      *  ANSI-C qsort()  *
  89      ********************/
  90
  91     lbStack[0] = (char *)base;
  92     ubStack[0] = (char *)base + (nmemb-1)*size;
  93     for (sp = 0; sp >= 0; sp--) {
  94         char *lb, *ub, *m;
  95         char *P, *i, *j;
  96
  97         lb = lbStack[sp];
  98         ub = ubStack[sp];
  99
 100         while (lb < ub) {
 101
 102             /* select pivot and exchange with 1st element */
 103             offset = (ub - lb) >> 1;
 104             P = lb + offset - offset % size;
 105             exchange (lb, P, size);
 106
 107             /* partition into two segments */
 108             i = lb + size;
 109             j = ub;
 110             while (1) {
 111                 while (i < j && compar(lb, i) > 0) i += size;
 112                 while (j >= i && compar(j, lb) > 0) j -= size;
 113                 if (i >= j) break;
 114                 exchange (i, j, size);
 115                 j -= size;
 116                 i += size;
 117             }
 118
 119             /* pivot belongs in A[j] */
 120             exchange (lb, j, size);
 121             m = j;
 122
 123             /* keep processing smallest segment, and stack largest */
 124             if (m - lb <= ub - m) {
 125                 if (m + size < ub) {
 126                     lbStack[sp] = m + size;
 127                     ubStack[sp++] = ub;
 128                 }
 129                 ub = m - size;
 130             } else {
 131                 if (m - size > lb) {
 132                     lbStack[sp] = lb;
 133                     ubStack[sp++] = m - size;
 134                 }
 135                 lb = m + size;
 136             }
 137         }
 138     }
 139 }
 140
 141 int pstrcmp(char *p1, char *p2);
 142
 143
 144
 145 /*********************************************************************/
 146 /*                                                                   */
 147 /*                                   HEAPSORT                           */
 148 /*                                                                   */
 149 /*                                                                   */
 150 /*********************************************************************/
 151 FILE * hsort_file ( FILE *fp,int n,int len, size_t size,int first_comp_field, int n_comp_fields,int (*compare)(const void *, const void*,int, int, size_t),void * (*copy)(void *,void*, size_t))
 152      {
 153      unsigned long i, ir, j, l;
 154      void *rra, *rrb, *ra_j, *ra_j_1;
 155      void *tp;
 156      long  start;
 157      int FE=1;
 158
 159
 160
 161      start=ftell(fp);
 162      rra   =vcalloc ( len, size);
 163      rrb   =vcalloc ( len, size);
 164      ra_j  =vcalloc ( len, size);
 165      ra_j_1=vcalloc ( len, size);
 166
 167      if ( n<2)return fp;
 168      l=(n >>1)+1;
 169      ir=n;
 170
 171      for (;;)
 172          {
 173          if ( l>FE)
 174             {
 175             l--;
 176             fseek( fp, start+(((l-1)*len)*size), SEEK_SET);
 177             fread( rra, size, len,fp); /*rra=ra[--l]*/
 178             }
 179          else
 180             {
 181             fseek( fp, start+((ir-1)*len*size), SEEK_SET);
 182             fread( rra, size, len,fp); /*rra=ra[ir]*/
 183
 184             fseek( fp, start, SEEK_SET);
 185             fread( rrb, size, len,fp); /*rrb=ra[0]*/
 186
 187             fseek( fp, start+((ir-1)*len*size), SEEK_SET);
 188             fwrite(rrb,size, len, fp); /*ra[ir]=rrb=ra[0]*/
 189
 190             if (--ir ==FE)
 191                {
 192                fseek ( fp,start, SEEK_SET);
 193                fwrite(rra,size, len, fp); /*ra[0]=rra*/
 194                break;
 195                }
 196             }
 197          i=l;
 198          j=l+l;
 199          while ( j<=ir)
 200                {
 201                fseek ( fp, start+((j-1)*len*size), SEEK_SET);
 202                fread (ra_j, size, len, fp);
 203
 204                if ( j<ir)
 205                   {
 206                   fseek ( fp, start+(((j-1)+1)*len*size), SEEK_SET);
 207                   fread (ra_j_1, size, len, fp);
 208                   }
 209
 210                if ( j<ir && compare( ra_j, ra_j_1,first_comp_field,n_comp_fields, size )<0)
 211                    {
 212                    SWAPP(ra_j, ra_j_1, tp);
 213                    j++;
 214                    }
 215
 216                if (compare(rra, ra_j, first_comp_field,n_comp_fields, size)<0)
 217                        {
 218                        fseek ( fp, start+((i-1)*len*size), SEEK_SET);
 219                        fwrite(ra_j,size, len, fp);
 220                        i=j;
 221                        j<<= 1;
 222                        }
 223                else
 224                        j=ir+1;
 225
 226                }
 227          fseek ( fp, start+((i-1)*len*size), SEEK_SET);
 228          fwrite(rra,size, len, fp);
 229          }
 230      vfree (  rra);
 231      vfree ( rrb);
 232
 233      vfree (  ra_j);
 234      vfree (  ra_j_1);
 235      return fp;
 236      }
 237 void ** hsort_array ( void **ra,int n,int len, size_t size,int first_comp_field, int n_comp_fields,int (*compare)(const void *, const void*,int, int, size_t),void * (*copy)(void *,void*, size_t))
 238      {
 239      unsigned long i, ir, j, l;
 240      void *rra;
 241      int FE=1;
 242
 243      if ( FE==1){ra--;}
 244      else
 245      {
 246      n--;
 247      }
 248
 249
 250      rra   =vcalloc ( len, size);
 251
 252
 253      if ( n<2)return ra;
 254      l=(n >>1)+1;
 255      ir=n;
 256
 257      for (;;)
 258          {
 259          if ( l>FE)
 260             {
 261             copy ( rra, ra[--l],len);
 262             }
 263          else
 264             {
 265             copy ( rra, ra[ir],len);
 266             copy ( ra[ir], ra[FE], len);
 267             if (--ir ==FE)
 268                {
 269                copy ( ra[FE],rra,len);
 270                break;
 271                }
 272             }
 273          i=l;
 274          j=l+l;
 275          while ( j<=ir)
 276                {
 277                if ( j<ir && compare( ra[j], ra[j+1],first_comp_field,n_comp_fields, size )<0)j++;
 278                if (compare(rra, ra[j], first_comp_field,n_comp_fields, size)<0)
 279                        {copy(ra[i], ra[j],len);i=j;j<<= 1;}
 280                else
 281                        j=ir+1;
 282                }
 283          copy( ra[i], rra,len);
 284          }
 285      vfree (rra);
 286      ra+=FE;
 287
 288      return ra;
 289      }
 290 /*********************************************************************/
 291 /*                                                                   */
 292 /*                         CEDRIC BSEARCH                            */
 293 /*                                                                   */
 294 /*                                                                   */
 295 /*********************************************************************/
 296 void * bsearch_file ( const void *key,int *p,int comp_first,int comp_len, FILE *fp,int len, int entry_len,size_t el_size, int (*compare)(const void *, const void*,int, int, size_t))
 297        {
 298        int upper, lower, c, i;
 299        static void *key2;
 300        long start;
 301        static long key2_size;
 302
 303
 304        start=ftell(fp);
 305
 306        upper=-1;
 307        lower=len;
 308        if ( key2==NULL){key2=vcalloc (entry_len, el_size);key2_size=entry_len* el_size;}
 309        else if (key2_size<  (entry_len* el_size)){vfree(key2);key2=vcalloc (entry_len, el_size);key2_size=entry_len* el_size;}
 310
 311        while ((lower-upper)>1)
 312              {
 313              i=(lower+upper) >> 1;
 314
 315              fseek ( fp,start+(i*el_size*entry_len), SEEK_SET);
 316              fread ( key2, el_size, entry_len,fp);
 317              c=compare(key2,key, comp_first, comp_len,el_size);
 318
 319              if      ( c==0){p[0]=i;return key2;}
 320              else if ( c< 0)upper=i;
 321              else if ( c> 0)lower=i;
 322              }
 323        return NULL;
 324        }
 325
 326 void * bsearch_array ( const void *key,int *p,int comp_first, int comp_len,void**list,int len, int entry_len,size_t el_size, int (*compare)(const void *, const void*,int, int, size_t))
 327        {
 328        int upper, lower, c, i;
 329        void *key2;
 330
 331        upper=-1;
 332        lower=len;
 333        while ((lower-upper)>1)
 334              {
 335              i=(lower+upper) >>1;
 336              key2=list[i];
 337              c=compare(key2,key, comp_first,comp_len,el_size);
 338
 339              if      ( c==0){p[0]=i;return key2;}
 340              else if ( c< 0)upper=i;
 341              else if ( c> 0)lower=i;
 342              }
 343        return NULL;
 344        }
 345
 346 /**********************************************************************/
 347 /*                                                                    */
 348 /*                         HSORT/BSEARCH WRAPPERS                             */
 349 /*                                                                    */
 350 /*                                                                    */
 351 /**********************************************************************/
 352 void **search_in_list_file ( void *key, int *p,int comp_len,FILE *fp, int len, size_t size, int entry_len)
 353       {
 354       static void **l;
 355
 356
 357       if ( l==NULL)l=vcalloc ( 1, sizeof (int*));
 358
 359       l[0]=bsearch_file (key,p,0,comp_len,fp,len,entry_len,size,hsort_cmp);
 360       if (l[0]==NULL)return NULL;
 361       else return l;
 362       }
 363 void **search_in_list_array ( void *key,int *p, int comp_len,void **L, int len, size_t size, int entry_len)
 364       {
 365       static void **l;
 366
 367       if ( l==NULL)l=vcalloc ( 1, sizeof (int*));
 368
 369       l[0]=bsearch_array (key,p,0,comp_len,L,len,entry_len,size,hsort_cmp);
 370       if (l[0]==NULL)return NULL;
 371       else return l;
 372       }
 373 void **hsort_list_array ( void **L, int len, size_t size, int entry_len, int first_comp_field, int n_comp_fields)
 374        {
 375          return hsort_array (L, len,entry_len, size,first_comp_field, n_comp_fields,hsort_cmp , hsort_cpy);
 376        }
 377 FILE  *hsort_list_file ( FILE*fp  , int len, size_t size, int entry_len, int first_comp_field, int n_comp_fields)
 378        {
 379
 380        return hsort_file (fp, len,entry_len, size,first_comp_field, n_comp_fields,hsort_cmp , hsort_cpy);
 381        }
 382
 383 int hsort_cmp ( const void *a, const void *b, int first, int clen, size_t size)
 384        {
 385        int*ax;
 386        int*bx;
 387        int p;
 388
 389        ax=(int*)a;
 390        bx=(int*)b;
 391        for ( p=first; p<clen+first; p++)
 392            {
 393            if ( ax[p]<bx[p])return -1;
 394            else if ( ax[p]==bx[p]);
 395            else return 1;
 396            }
 397        return 0;
 398        }
 399 void *hsort_cpy(void*to, void *from, size_t size)
 400        {
 401
 402        int *ax;
 403        int *bx;
 404        int p;
 405        ax=(int*)to;
 406        bx=(int*)from;
 407        for (p=0; p<(int)size; p++)
 408            ax[p]=bx[p];
 409
 410
 411        return to;
 412
 413        }
 414
 415
 416 void test_hsort_list_array()
 417       {
 418       int **array;
 419       int a;
 420       int n=100;
 421
 422       array=declare_int(n, 3);
 423       for ( a=0; a<n; a++)array[a][0]=a;
 424
 425       hsort_list_array( (void**)array,n, sizeof (int), 3, 0, 1);
 426       for ( a=0; a<n; a++)fprintf ( stderr, "\n%d %d", array[a][0],a);
 427       myexit(EXIT_FAILURE);
 428       }
 429
 430
 431 /*********************************************************************/
 432 /*                                                                   */
 433 /*                         B_SEARCH_FILE FUNCTIONS                   */
 434 /*                                                                   */
 435 /*                                                                   */
 436 /*********************************************************************/
 437
 438 /*********************************************************************/
 439 /*                                                                   */
 440 /*                         SORT/COMPARE/SEARCH FUNCTIONS             */
 441 /*                                                                   */
 442 /*                                                                   */
 443 /*********************************************************************/
 444 static int sort_field;
 445 int **search_in_list_int ( int *key, int k_len, int **list, int ne)
 446         {
 447         int **l;
 448         sort_field=k_len;
 449         l=bsearch (&key,list, ne, sizeof(int**),(int(*)(const void*,const void*))(cmp_list_int));
 450         return l;
 451         }
 452 void sort_float ( float **V,int N_F, int F, int left, int right)
 453         {
 454         sort_field=F;
 455         qsort ( V, right+1, sizeof(int**),(int(*)(const void*,const void*))(cmp_int));
 456         }
 457 int cmp_float ( const float **a, const float **b)
 458         {
 459         if ( a[0][sort_field]< b[0][sort_field])return-1;
 460         else if ( a[0][sort_field]==b[0][sort_field])return 0;
 461         else return 1;
 462         }
 463
 464 void sort_int_1D ( int *L, int n)
 465         {
 466         int **array;
 467         int a;
 468
 469         array=declare_int ( n, 1);
 470         for ( a=0; a< n; a++)
 471                 array[a][0]=L[a];
 472         sort_int ( array, 1, 0, 0, n-1);
 473         for ( a=0; a< n; a++)
 474                 L[a]=array[a][0];
 475         free_int ( array, n);
 476         }
 477
 478 char** sort_string_array (char **V, int n)
 479 {
 480
 481
 482   qsort ( V,n,sizeof (char*),(int(*)(const void*,const void*))(pstrcmp));
 483   return V;
 484
 485
 486 }
 487
 488 int pstrcmp(char *p1, char *p2)
 489 {
 490   return strcmp(*(char **)p1, *(char **)p2);
 491 }
 492 void sort_int ( int **V,int N_F, int F, int left, int right)
 493         {
 494           if (!V)return;
 495         sort_field=F;
 496         qsort ( V, (right-left)+1, sizeof(int**),(int(*)(const void*,const void*))(cmp_int));
 497         }
 498 void sort_list_int ( int **V,int N_F, int F, int left, int right)
 499         {
 500         if (!V)return;
 501         sort_field=F;
 502         qsort ( V, (right-left)+1, sizeof(int**),(int(*)(const void*,const void*))(cmp_list_int));
 503         }
 504 static int *order;
 505 void sort_list_int2 ( int **V,int *list,int N_F, int left, int right)
 506         {
 507           // just like sort_int_list, but uses list to to order the comparison of the keys
 508           if (!V)return;
 509         order=list;
 510
 511         qsort ( V, (right-left)+1, sizeof(int**),(int(*)(const void*,const void*))(cmp_list_int2));
 512         }
 513 int cmp_list_int2 (const int**a, const int**b)
 514         {
 515           int p=0;;
 516           int c;
 517
 518
 519
 520           while ((c=order[p])!=-1)
 521           {
 522             if ( a[0][c]>b[0][c])return   1;
 523             else if ( a[0][c]<b[0][c])return  -1;
 524             p++;
 525           }
 526         return 0;
 527         }
 528 void sort_int_inv ( int **V,int N_F, int F, int left, int right)
 529         {
 530         int a,b;
 531         int **list;
 532         if (!V)return;
 533         sort_field=F;
 534         qsort ( V, (right-left)+1, sizeof(int**),(int(*)(const void*,const void*))(cmp_int));
 535
 536         list=declare_int ((right-left)+1, N_F);
 537         for ( a=left; a< (right-left)+1; a++)
 538                 {
 539                 for ( b=0; b< N_F; b++)
 540                         {
 541                         list[a-left][b]=V[a][b];
 542                         }
 543                 }
 544         for ( a=left; a< (right-left)+1; a++)
 545                 {
 546                 for ( b=0; b< N_F; b++)
 547                         V[a][b]=list[(right-left)-a][b];
 548                 }
 549         free_int (list, -1);
 550         }
 551
 552 void sort_list_int_inv ( int **V,int N_F, int F, int left, int right)
 553         {
 554         int a,b;
 555         int **list;
 556         if (!V)return;
 557         sort_field=F;
 558         qsort ( V, (right-left)+1, sizeof(int**),(int(*)(const void*,const void*))(cmp_list_int));
 559
 560         list=declare_int ((right-left)+1, N_F);
 561         for ( a=left; a< (right-left)+1; a++)
 562                 {
 563                 for ( b=0; b< N_F; b++)
 564                         {
 565                         list[a-left][b]=V[a][b];
 566                         }
 567                 }
 568         for ( a=left; a< (right-left)+1; a++)
 569                 {
 570                 for ( b=0; b< N_F; b++)
 571                         V[a][b]=list[(right-left)-a][b];
 572                 }
 573         free_int (list, -1);
 574         }
 575
 576
 577
 578 int cmp_int ( const int**a, const int**b)
 579         {
 580         if ( a[0][sort_field]< b[0][sort_field])return-1;
 581         else if ( a[0][sort_field]==b[0][sort_field])return 0;
 582         else return 1;
 583         }
 584 int cmp_list_int (const int**a, const int**b)
 585         {
 586         int c;
 587         int undef=0;
 588         int ret;
 589
 590         for ( c=0; c<=sort_field; c++)
 591                 {
 592                 if ( a[0][c]==UNDEFINED|| b[0][c]==UNDEFINED)ret=0;
 593                 else if ( a[0][c]>b[0][c])return   1;
 594                 else if ( a[0][c]<b[0][c])return  -1;
 595                 }
 596         if (undef==sort_field)
 597                 {
 598                   if (a[0][0]==b[0][0])return 0;
 599                 }
 600         return 0;
 601         }
 602
 603
 604 int name_is_in_list ( char *name, char **name_list, int n_name, int len)
 605         {
 606         int a;
 607         int pos=-1;
 608         /*Note: RETURNS THE Offset of the LAST Occurence of name in name_list*/
 609
 610         if ( name_list==NULL || name ==NULL)return -1;
 611
 612
 613         for ( a=0; a< n_name; a++)
 614           {
 615             if ( name_list[a]==NULL);
 616             else if ( len!=-1)
 617               {
 618                 if (strncmp ( name, name_list[a], len)==0)pos=a;
 619               }
 620             else if ( strm ( name, name_list[a]))pos=a;
 621           }
 622         return pos;
 623         }
 624 char * check_list_for_dup ( char **list, int ne)
 625         {
 626         int a, b;
 627
 628         for ( a=0; a< ne-1; a++)
 629             for ( b=a+1; b< ne; b++)if (strm ( list[a], list[b]))return list[a];
 630         return NULL;
 631         }
 632 FILE *get_number_list_in_file ( FILE *fp, int *list, int *n, int *max_len)
 633         {
 634
 635         int c;
 636
 637         while ( isspace((c=fgetc (fp))));
 638         ungetc(c, fp);
 639         while ( c!='\n')
 640                 {
 641                 while ( isspace((c=fgetc (fp))) && c!='\n');
 642
 643                 if ( c!='\n')
 644                         {
 645                         ungetc(c, fp);
 646                         if ( n[0]>=max_len[0])
 647                                 list=vrealloc ( list, (n[0]+100)*sizeof (int));
 648                                 max_len[0]=(n[0]+100);
 649
 650                         fscanf ( fp, "%d",&list[n[0]++]);
 651                         }
 652                 }
 653         return fp;
 654         }
 655 /*********************************************************************/
 656 /*                                                                   */
 657 /*                         QUANTILE                                  */
 658 /*                                                                   */
 659 /*                                                                   */
 660 /*********************************************************************/
 661 int quantile (int argc, char *argv[])
 662 {
 663   FILE *fp;
 664   int a,n,c, t;
 665   int **list;
 666   char ***  string_list;
 667   char *name, s1[1000], s2[1000];
 668
 669
 670   if ( argc<2)
 671     {
 672
 673       fprintf (stderr, "\nquantile <fname> <quant: 0.00-1.00> [<top | bottom>]");
 674       fprintf (stderr, "\nSplits your data in two according to the quantile");
 675       fprintf (stderr, "\nReturns the top quantile or the bottom quantile");
 676       fprintf (stderr, "\nData must be in <fname> with two fields/line: Field1=index, Field2=value\n");
 677       fprintf (stderr, "\n1 27\n2 46\n3 5\n...\n");
 678       fprintf (stderr, "\nValue can either be integer or float");
 679
 680
 681       myexit (EXIT_FAILURE);
 682     }
 683
 684   if (strm (argv[1], "stdin"))
 685     {
 686       name=vtmpnam(NULL);
 687       fp=vfopen (name, "w");
 688       while ( (c=fgetc(stdin))!=EOF)
 689         {
 690           fprintf ( fp, "%c", c);
 691         }
 692       vfclose (fp);
 693     }
 694   else
 695     name=argv[1];
 696
 697
 698
 699
 700   n=count_n_line_in_file (name);
 701   list=declare_int (n, 2);
 702   string_list=declare_arrayN(3,sizeof (char), n, 2, 10);
 703
 704   fp=vfopen (name, "r");
 705   n=0;
 706   while ( (c=fgetc (fp))!=EOF)
 707     {
 708       ungetc(c,fp);
 709       fscanf ( fp, "%s %s\n", s1, s2);
 710       list[n][0]=(int)(atof(s1)*1000);
 711       list[n][1]=(int)(atof(s2)*1000);
 712       list[n][2]=n;
 713
 714       sprintf (string_list[n][0],"%s",s1);
 715       sprintf (string_list[n][1],"%s",s2);
 716
 717
 718       n++;
 719     }
 720   sort_int_inv ( list,3, 1, 0, n-1);
 721   t=quantile_rank ( list,1,n, atof (argv[2]));
 722   if ( argc!=4 || (argc==4 && strm (argv[3], "bottom")))
 723     {
 724       for (a=t; a<n; a++)
 725         fprintf ( stdout, "%s %s\n", string_list[list[a][2]][0], string_list[list[a][2]][1]);
 726     }
 727   else
 728     {
 729       for ( a=0; a<t; a++)
 730         fprintf ( stdout, "%s %s\n", string_list[list[a][2]][0], string_list[list[a][2]][1]);
 731     }
 732
 733   fprintf (stderr, "\nQuantile %.2f T= %d Out of N= %d entries\n", atof (argv[2]), t, n),
 734   free_int (list, -1);
 735   return n;
 736
 737 }
 738
 739 int quantile_rank (int **list,int field, int n, float p)
 740 {
 741   int nr;
 742
 743
 744   if ( p==1) nr=0;
 745   else if ( p==0) nr=n;
 746   else
 747     {
 748       int a, b,j, *l;
 749       double g, q, np, i_part;
 750       l=vcalloc ( n, sizeof (int));
 751       for (a=n-1, b=0; b<n; a--, b++)
 752         l[b]=list[a][field];
 753
 754       np=(double)n*(double)p;
 755       g=modf (np, &i_part);
 756       j=(int)i_part;
 757       j--;
 758
 759       q=(float)l[j]+g*((float)l[j+1]-(float)l[j]);
 760
 761       nr=0;
 762       while (nr<n && list[nr][field]>=q)nr++;
 763       vfree(l);
 764     }
 765   return nr;
 766 }
 767 /*********************************************************************/
 768 /*                                                                   */
 769 /*                         DUPLICATION                               */
 770 /*                                                                   */
 771 /*                                                                   */
 772 /*********************************************************************/
 773 short* ga_memcpy_short ( short *array1, short *array2, int n)
 774         {
 775         int a;
 776
 777         for ( a=0; a< n; a++)
 778                 array2[a]=array1[a];
 779
 780         return array2;
 781         }
 782 int * ga_memcpy_int ( int *array1, int *array2, int n)
 783         {
 784         int a;
 785
 786         for ( a=0; a< n; a++)
 787                 array2[a]=array1[a];
 788
 789         return array2;
 790         }
 791
 792 float* ga_memcpy_float ( float *array1, float *array2, int n)
 793         {
 794         int a;
 795
 796         for ( a=0; a< n; a++)
 797                 array2[a]=array1[a];
 798
 799         return array2;
 800         }
 801 double*  ga_memcpy_double (double *array1, double*array2, int n)
 802         {
 803         int a;
 804
 805         for ( a=0; a< n; a++)
 806                 array2[a]=array1[a];
 807
 808         return array2;
 809         }
 810
 811
 812
 813 /*recycle: get the bottom pointer on the top of the heap*/
 814
 815 void ** recycle (void **A, int l, int cycle)
 816 {
 817   void **B;
 818   int a,b,c;
 819   B=vcalloc (l, sizeof (void*));
 820
 821   for ( c=0; c< cycle; c++)
 822     {
 823       for ( a=1, b=0; a<l; a++, b++) B[b]=A[a];
 824       B[l-1]=A[0];
 825       for ( a=0; a<l; a++)A[a]=B[a];
 826     }
 827   vfree (B);
 828   return A;
 829 }
 830
 831 /* Old READ/WRITE ARRAY SIZE*/
 832 /*
 833 #define WRITE_SIZE(type,function)\
 834 void function ( int x, type *array, int os)\
 835      {\
 836      int a,l;\
 837      char buf[SIZE_OF_INT+1];\
 838      array+=os*SIZE_OF_INT;\
 839      for ( a=0;a<SIZE_OF_INT; a++)array[a]=0;\
 840      sprintf ( buf, "%d", x);\
 841      l=strlen (buf);\
 842      array+=SIZE_OF_INT-l;\
 843      for (a=0; a<l; a++)array[a]=(type)buf[a];\
 844      }
 845 WRITE_SIZE(short,write_size_short)
 846 WRITE_SIZE(char,write_size_char)
 847 WRITE_SIZE(int,write_size_int)
 848 WRITE_SIZE(float,write_size_float)
 849 WRITE_SIZE(double,write_size_double)
 850
 851 #define READ_ARRAY_SIZE(type, function)\
 852 int function (type *array, int os)\
 853     {\
 854     int a, b;\
 855     char buf[SIZE_OF_INT+1];\
 856     a=b=0;\
 857     array+=os*SIZE_OF_INT;\
 858     while ( a!=SIZE_OF_INT && array[a]==0)a++;\
 859     while ( a!=SIZE_OF_INT)buf[b++]=(char)array[a++];\
 860     buf[b]='\0';\
 861     return atoi(buf);\
 862     }
 863 READ_ARRAY_SIZE(short,read_size_short)
 864 READ_ARRAY_SIZE(char,read_size_char)
 865 READ_ARRAY_SIZE(int,read_size_int)
 866 READ_ARRAY_SIZE(float,read_size_float)
 867 READ_ARRAY_SIZE(double,read_size_double)
 868 */
 869 /*********************************************************************/
 870 /*                                                                   */
 871 /*                          DUPLICATION                              */
 872 /*                                                                   */
 873 /*                                                                   */
 874 /*********************************************************************/
 875
 876 #define SET_NUMBERS(type,function)\
 877 type * function(type *list,int n, ...)\
 878      {\
 879      type *buf;\
 880      int  *index,i;\
 881      va_list ap;\
 882      int max;\
 883 \
 884      va_start(ap, n);\
 885      buf=vcalloc(n, sizeof(type));\
 886      index=vcalloc(n, sizeof(int));\
 887      max=0;\
 888      for ( i=0; i< n; i++)\
 889          {\
 890          buf[i]  =va_arg(ap,type);\
 891          index[i]=va_arg(ap, int);\
 892          if (index[i]>max)max=index[i];\
 893          }\
 894      va_end(ap);\
 895      if (list==NULL)list=vcalloc ( max+1, sizeof (type));\
 896      for ( i=0; i<n; i++)list[index[i]]=buf[i];\
 897      vfree(buf);\
 898      vfree(index);\
 899      return list;\
 900      }
 901      /*SET_NUMBERS(short ,set_short)*/
 902      /*SET_NUMBERS(char  ,set_char)*/
 903 SET_NUMBERS(int   ,set_int)
 904      /*SET_NUMBERS(float ,set_float)*/
 905 SET_NUMBERS(double,set_double)
 906
 907 short ** duplicate_short ( short **array , int len, int field)
 908     {
 909     return copy_short (array ,declare_short ( len, field),  len, field);
 910     }
 911 int ** duplicate_int ( int **array , int len, int field)
 912     {
 913     return copy_int (array ,declare_int ( len, field),  len, field);
 914     }
 915 char ** duplicate_char ( char **array , int len, int field)
 916     {
 917     return copy_char (array ,declare_char ( len, field),  len, field);
 918     }
 919 char * duplicate_string ( char *string)
 920     {
 921       int l;
 922       char *buf=NULL;
 923
 924       l=strlen (string);
 925
 926       if ( !l);
 927       else
 928         {
 929           buf=vcalloc ( l+1, sizeof(char));
 930           sprintf ( buf, "%s", string);
 931         }
 932       return buf;
 933     }
 934 float ** duplicate_float ( float **array , int len, int field)
 935     {
 936     return copy_float (array ,declare_float ( len, field),  len, field);
 937     }
 938 double ** duplicate_double ( double **array , int len, int field)
 939     {
 940     return copy_double (array ,declare_double ( len, field),  len, field);
 941     }
 942
 943
 944
 945 /*********************************************************************/
 946 /*                                                                   */
 947 /*                           COPY OF 2D ARRAY                        */
 948 /*                                                                   */
 949 /*                                                                   */
 950 /*********************************************************************/
 951 short ** copy_short( short **array1, short  **array2, int len, int number_field)
 952     {
 953     int a;
 954
 955     if ( len==-1)len=read_array_size (array1,sizeof (short*));
 956     if ( number_field==-1)number_field=read_array_size (array1[0],sizeof (short));
 957     if ( array2)free_short ( array2, -1);
 958     array2=declare_short ( len, number_field);
 959
 960     for ( a=0; a< len; a++)
 961         ga_memcpy_short( array1[a],array2[a],number_field);
 962
 963     return array2;
 964     }
 965 char ** copy_char ( char **array1, char **array2, int len, int number_field)
 966     {
 967     int a;
 968
 969     if ( array1==NULL)return NULL;
 970     if ( len==-1)len=read_size_char (array1,sizeof (char*));
 971     if ( number_field==-1)
 972         {
 973         number_field=read_size_char (array1[0],sizeof(char));
 974         for ( a=0; a< len; a++)
 975             number_field=MAX(number_field, strlen ( array1[a]))+1;
 976         }
 977
 978     if ( array2)free_char (array2, -1);
 979     array2=declare_char(len, number_field);
 980
 981     for ( a=0; a< len; a++)
 982       sprintf ( array2[a], "%s", array1[a]);
 983
 984     return array2;
 985     }
 986 int ** copy_int ( int **array1, int **array2, int len, int number_field)
 987     {
 988     int a;
 989
 990     if ( array1==NULL) return NULL;
 991     if ( len==-1)len=read_size_int (array1, sizeof (int*));
 992     if ( number_field==-1)number_field=read_size_int (array1[0],sizeof (int));
 993
 994
 995
 996     if (array2)free_int (array2, -1);
 997     array2=declare_int ( len, number_field);
 998
 999     for ( a=0; a< len; a++)
1000         ga_memcpy_int( array1[a],array2[a],number_field);
1001
1002     return array2;
1003     }
1004
1005 float ** copy_float ( float **array1, float **array2, int len, int number_field)
1006     {
1007     int a;
1008
1009     if ( array1==NULL) return NULL;
1010     if ( len==-1)len=read_size_float (array1,sizeof (float*));
1011     if ( number_field==-1)number_field=read_size_float (array1[0],sizeof (float));
1012
1013     if ( array2)free_float (array2, -1);
1014     array2=declare_float ( len, number_field);
1015
1016     for ( a=0; a< len; a++)
1017         ga_memcpy_float( array1[a],array2[a],number_field);
1018     return array2;
1019     }
1020 double ** copy_double (double **array1, double **array2, int len, int number_field)
1021     {
1022     int a;
1023
1024     if ( array1==NULL) return NULL;
1025     if ( len==-1)len=read_size_double (array1,sizeof (double*));
1026     if ( number_field==-1)number_field=read_size_double (array1[0], sizeof (double));
1027     if ( array2)free_double (array2, -1);
1028     array2=declare_double ( len, number_field);
1029
1030     for ( a=0; a< len; a++)
1031         ga_memcpy_double( array1[a],array2[a],number_field);
1032     return array2;
1033     }
1034 /*********************************************************************/
1035 /*                                                                   */
1036 /*                        CONCATENATION                              */
1037 /*                                                                   */
1038 /*                                                                   */
1039 /*********************************************************************/
1040
1041
1042
1043
1044
1045 Alignment ** cat_aln_list ( Alignment **list_to_cat,int first, int end, Alignment **rec_list)
1046     {
1047     int rec_list_start;
1048     int a, b;
1049
1050     if ( list_to_cat==NULL)return rec_list;
1051     else
1052        {
1053        rec_list_start=(rec_list[-1])->nseq;
1054        rec_list=realloc_aln_array ( rec_list, end-first);
1055        for ( a=first, b=rec_list_start; a<end; a++, b++)copy_aln (list_to_cat[a], rec_list[b]);
1056        free_aln_array ( list_to_cat);
1057        return rec_list;
1058        }
1059     }
1060
1061 /*********************************************************************/
1062 /*                                                                   */
1063 /*                         NUMBER ARRAY ANALYSE                      */
1064 /*                                                                   */
1065 /*                                                                   */
1066 /*********************************************************************/
1067
1068 FILE * output_array_int (int  **array, int len, int nf ,FILE *fp)
1069      {
1070      int a, b;
1071      for ( a=0; a<len; a++)
1072          {fprintf ( fp, "\n");
1073             for (b=0; b<nf; b++)fprintf (fp, "%5d ", array[a][b]);
1074
1075          }
1076      fprintf ( fp, "\n");
1077      return fp;
1078      }
1079 #define RETURN_MAX_COOR(type,wf,rf,function,comparison,undef)\
1080 type function ( type ** array, int len_array, int field, int *coor)\
1081     {\
1082     type max;\
1083     int a;\
1084 \
1085     if (array==NULL || len_array==0)return 0;\
1086     else\
1087         {\
1088         len_array=rf(array,sizeof (type*));\
1089         max=array[0][field];\
1090         coor[0]=0;\
1091         for ( a=1; a< len_array; a++)\
1092               if ( max==undef)max=array[a][field];\
1093               else if ( array[a][field]!=undef)\
1094                    if (array[a][field] comparison max)\
1095                        {max=array[a][field];\
1096                         coor[0]=a;\
1097                         }\
1098          }\
1099     return max;\
1100     }
1101
1102 RETURN_MAX_COOR(short,write_size_short,read_size_short,return_max_coor_short,>, UNDEFINED_SHORT)
1103 RETURN_MAX_COOR(char,write_size_char,read_size_char,return_max_coor_char,>, UNDEFINED_CHAR)
1104 RETURN_MAX_COOR(int,write_size_int,read_size_int,return_max_coor_int,>, UNDEFINED_INT)
1105 RETURN_MAX_COOR(float,write_size_float,read_size_float,return_max_coor_float,>, UNDEFINED_FLOAT)
1106 RETURN_MAX_COOR(double,write_size_double,read_size_double,return_max_coor_double,>, UNDEFINED_DOUBLE)
1107 RETURN_MAX_COOR(short,write_size_short,read_size_short,return_min_coor_short,<, UNDEFINED_SHORT)
1108 RETURN_MAX_COOR(char,write_size_char,read_size_char,return_min_coor_char,<, UNDEFINED_CHAR)
1109 RETURN_MAX_COOR(int,write_size_int,read_size_int,return_min_coor_int,<, UNDEFINED_INT)
1110 RETURN_MAX_COOR(float,write_size_float,read_size_float,return_min_coor_float,<, UNDEFINED_FLOAT)
1111 RETURN_MAX_COOR(double,write_size_double,read_size_double,return_min_coor_double,<, UNDEFINED_DOUBLE)
1112 #define RETURN_MAX(type,wf,rf,function,comparison,undef)\
1113 type function ( type ** array, int len_array, int field)\
1114     {\
1115     type max;\
1116     int a;\
1117 \
1118     if (array==NULL || len_array==0)return 0;\
1119     else\
1120         {\
1121         if (len_array==-1)len_array=rf(array,sizeof (type*));\
1122         max=array[0][field];\
1123         for ( a=1; a< len_array; a++)\
1124             if ( max==undef)max=array[a][field];\
1125             else if ( array[a][field]!=undef)max=( array[a][field] comparison max)?array[a][field]:max;\
1126         }\
1127     return (max==undef)?0:max;\
1128     }
1129
1130 RETURN_MAX(short,write_size_short,read_size_short,return_max_short,>,UNDEFINED_SHORT)
1131 RETURN_MAX(char,write_size_char,read_size_char,return_max_char,>,UNDEFINED_CHAR)
1132 RETURN_MAX(int,write_size_int,read_size_int,return_max_int,>,UNDEFINED_INT)
1133 RETURN_MAX(float,write_size_float,read_size_float,return_max_float,>,UNDEFINED_FLOAT)
1134 RETURN_MAX(double,write_size_double,read_size_double,return_max_double,>,UNDEFINED_DOUBLE)
1135 RETURN_MAX(short,write_size_short,read_size_short,return_min_short,<,UNDEFINED_SHORT)
1136 RETURN_MAX(char,write_size_char,read_size_char,return_min_char,<,UNDEFINED_CHAR)
1137 RETURN_MAX(int,write_size_int,read_size_int,return_min_int,<,UNDEFINED_INT)
1138 RETURN_MAX(float,write_size_float,read_size_float,return_min_float,<,UNDEFINED_FLOAT)
1139 RETURN_MAX(double,write_size_double,read_size_double,return_min_double,<,UNDEFINED_DOUBLE)
1140
1141
1142
1143 #define RETURN_2DMAX(type,wf,rf,function,comparison,undef)\
1144 type function ( type ** array, int start, int len_array, int first_field, int number_field)\
1145     {\
1146     type max;\
1147     int a,b;\
1148     if (array==NULL || len_array==0 || first_field<0 || number_field==0)return 0;\
1149     else\
1150          {max=array[start][first_field];\
1151           for ( a=start; a< start+len_array; a++)\
1152               for (b=first_field; b< first_field+number_field; b++)\
1153                      if (array[a][b]!=undef)max=( array[a][b] comparison max)?array[a][b]:max;\
1154          }\
1155     return max;\
1156     }
1157 RETURN_2DMAX(short,write_size_short,read_size_short,return_2Dmax_short,>, UNDEFINED_SHORT)
1158 RETURN_2DMAX(char,write_size_char,read_size_char,return_2Dmax_char,>,UNDEFINED_CHAR)
1159 RETURN_2DMAX(int,write_size_int,read_size_int,return_2Dmax_int,>,UNDEFINED_INT)
1160 RETURN_2DMAX(float,write_size_float,read_size_float,return_2Dmax_float,>,UNDEFINED_FLOAT)
1161 RETURN_2DMAX(double,write_size_double,read_size_double,return_2Dmax_double,>,UNDEFINED_DOUBLE)
1162 RETURN_2DMAX(short,write_size_short,read_size_short,return_2Dmin_short,<,UNDEFINED_SHORT)
1163 RETURN_2DMAX(char,write_size_char,read_size_char,return_2Dmin_char,<,UNDEFINED_CHAR)
1164 RETURN_2DMAX(int,write_size_int,read_size_int,return_2Dmin_int,<,UNDEFINED_INT)
1165 RETURN_2DMAX(float,write_size_float,read_size_float,return_2Dmin_float,<,UNDEFINED_FLOAT)
1166 RETURN_2DMAX(double,write_size_double,read_size_double,return_2Dmin_double,<,UNDEFINED_DOUBLE)
1167
1168 #define RETURN_2DMAX_COOR(type,wf,rf,function,compare,undef)\
1169 type function ( type **array, int start1 , int end1, int start2, int end2,int *i, int *j)\
1170     {\
1171     int a, b;\
1172     double max=undef;\
1173     if ( start1==-1)start1=0;\
1174     if ( start2==-1)start2=0;\
1175     if ( end1==-1)end1=rf(array,sizeof (type*));\
1176     if ( end2==-1)end2=rf(array[0],sizeof (type));\
1177     if ( array==NULL || (end1-start1)==0 || (end1-start1)>rf ( array,sizeof (type*)) || (end2-start2)==0)\
1178         {\
1179         return 0;\
1180         i[0]=0;\
1181         j[0]=0;\
1182         }\
1183     i[0]=0;\
1184     j[0]=0;\
1185     for ( a=start1; a<end1; a++)\
1186         for ( b=start2; b<end2; b++)\
1187             {\
1188             if ( max==undef && array[a][b]!=undef)max=array[a][b];\
1189             else if ( array[a][b]!=undef && (array[a][b] compare max))\
1190                {\
1191                max=array[a][b];\
1192                i[0]=a;\
1193                j[0]=b;\
1194                }\
1195             }\
1196     return (type)max;\
1197     }
1198 RETURN_2DMAX_COOR(short,write_size_short,read_size_short,return_2Dmax_coor_short,>,UNDEFINED_SHORT)
1199 RETURN_2DMAX_COOR(char,write_size_char,read_size_char,return_2Dmax_coor_char,>,UNDEFINED_CHAR)
1200 RETURN_2DMAX_COOR(int,write_size_int,read_size_int,return_2Dmax_coor_int,>,UNDEFINED_INT)
1201 RETURN_2DMAX_COOR(float,write_size_float,read_size_float,return_2Dmax_coor_float,>,UNDEFINED_FLOAT)
1202 RETURN_2DMAX_COOR(double,write_size_double,read_size_double,return_2Dmax_coor_double,>,UNDEFINED_DOUBLE)
1203 RETURN_2DMAX_COOR(short,write_size_short,read_size_short,return_2Dmin_coor_short,<,UNDEFINED_SHORT)
1204 RETURN_2DMAX_COOR(char,write_size_char,read_size_char,return_2Dmin_coor_char,<,UNDEFINED_CHAR)
1205 RETURN_2DMAX_COOR(int,write_size_int,read_size_int,return_2Dmin_coor_int,<,UNDEFINED_INT)
1206 RETURN_2DMAX_COOR(float,write_size_float,read_size_float,return_2Dmin_coor_float,<,UNDEFINED_FLOAT)
1207 RETURN_2DMAX_COOR(double,write_size_double,read_size_double,return_2Dmin_coor_double,<,UNDEFINED_DOUBLE)
1208
1209 #define RETURN_WMEAN(type,wf,rf,function,sum_function,undef)\
1210 double function ( type **array, int len, int wfield,int sfield)\
1211     {\
1212     double b;\
1213     int a, c;\
1214     if ( len==0 ||array==NULL || len>rf ( array,sizeof (type*)))return 0;\
1215     else\
1216          {\
1217          if ( len==-1)len=rf(array,sizeof (type*));\
1218          for ( b=0, c=0,a=0; a< len; a++)\
1219              {\
1220              if (array[a][sfield]!=undef && array[a][wfield]!=undef )\
1221                 {\
1222                 b+=array[a][sfield];\
1223                 c+=array[a][wfield];\
1224                 }\
1225              }\
1226          }\
1227     return (c==0)?0:(b/c);\
1228     }
1229 RETURN_WMEAN(short,write_size_short,read_size_short,return_wmean_short, return_sum_short,UNDEFINED_SHORT)
1230 RETURN_WMEAN(char,write_size_char,read_size_char, return_wmean_char,return_sum_char,UNDEFINED_CHAR)
1231 RETURN_WMEAN(int,write_size_int,read_size_int,return_wmean_int,return_sum_int,UNDEFINED_INT)
1232 RETURN_WMEAN(float,write_size_float,read_size_float,return_wmean_float,return_sum_float,UNDEFINED_FLOAT)
1233 RETURN_WMEAN(double,write_size_double,read_size_double,return_wmean_double,return_sum_double,UNDEFINED_DOUBLE)
1234
1235
1236 #define RETURN_MEAN(type,wf,rf,function,sum_function,undef)\
1237 double function ( type **array, int len, int field)\
1238     {\
1239     double b;\
1240     int a, c;\
1241     if ( len==0 ||array==NULL || len>rf ( array,sizeof(type*)))return 0;\
1242     else\
1243          {\
1244          for ( b=0, c=0,a=0; a< len; a++)\
1245              {\
1246              if (array[a][field]!=undef)\
1247                 {\
1248                 b+=array[a][field];\
1249                 c++;\
1250                 }\
1251              }\
1252          }\
1253     return (c==0)?0:(b/c);\
1254     }
1255 RETURN_MEAN(short,write_size_short,read_size_short,return_mean_short, return_sum_short,UNDEFINED_SHORT)
1256 RETURN_MEAN(char,write_size_char,read_size_char, return_mean_char,return_sum_char,UNDEFINED_CHAR)
1257 RETURN_MEAN(int,write_size_int,read_size_int,return_mean_int,return_sum_int,UNDEFINED_INT)
1258 RETURN_MEAN(float,write_size_float,read_size_float,return_mean_float,return_sum_float,UNDEFINED_FLOAT)
1259 RETURN_MEAN(double,write_size_double,read_size_double,return_mean_double,return_sum_double,UNDEFINED_DOUBLE)
1260
1261 #define RETURN_SUM(type,wf,rf,function,undef)\
1262 type function(type **array, int len, int field)\
1263 {\
1264  int a;\
1265  type b=0;\
1266  if ( len==0 ||array==NULL)return 0;\
1267  else\
1268      {\
1269      if ( len==-1)len=rf ( array,sizeof (type*));\
1270      for ( a=0; a< len; a++)\
1271           if ( array[a][field]!=undef)b+=array[a][field];\
1272      }\
1273   return b;\
1274   }
1275 RETURN_SUM(short,write_size_short,read_size_short, return_sum_short,UNDEFINED_SHORT)
1276 RETURN_SUM(char,write_size_char,read_size_char,return_sum_char,UNDEFINED_CHAR)
1277 RETURN_SUM(int,write_size_int,read_size_int,return_sum_int,UNDEFINED_INT)
1278 RETURN_SUM(float,write_size_float,read_size_float,return_sum_float,UNDEFINED_FLOAT)
1279 RETURN_SUM(double,write_size_double,read_size_double,return_sum_double,UNDEFINED_DOUBLE)
1280
1281 #define RETURN_SD(type,wf,rf,function,undef)\
1282   type function ( type **array, int len, int field,type mean)   \
1283     {\
1284     int a;\
1285     double c=0;\
1286     if ( len==0 ||array==NULL || len>rf ( array,sizeof(type*)))return 0;\
1287     else\
1288         {\
1289         for ( a=0; a< len; a++)\
1290             {\
1291             if ((array[a][field]!=undef) && (mean-array[a][field])!=0)\
1292              c+=((double)mean-array[a][field])*((double)mean-array[a][field]);\
1293             }\
1294         c=sqrt(c)/(double)len;\
1295         return (type)MAX(c,1);\
1296         }\
1297     }
1298 RETURN_SD(short,write_size_short,read_size_short, return_sd_short,UNDEFINED_SHORT)
1299 RETURN_SD(char,write_size_char,read_size_char,return_sd_char,UNDEFINED_CHAR)
1300 RETURN_SD(int,write_size_int,read_size_int,return_sd_int,UNDEFINED_INT)
1301 RETURN_SD(float,write_size_float,read_size_float,return_sd_float,UNDEFINED_FLOAT)
1302 RETURN_SD(double,write_size_double,read_size_double,return_sd_double,UNDEFINED_DOUBLE)
1303 double return_z_score( double x,double sum, double sum2, double n)
1304     {
1305     double sd;
1306     double avg;
1307     double z;
1308
1309
1310     sd=(n==0)?0:sqrt(sum2*n -sum*sum)/n;
1311     avg=(n==0)?0:(sum/n);
1312     z=(sd==0)?0:(x-avg)/sd;
1313     return z;
1314     }
1315
1316 double* return_r (double **list, int n)
1317 {
1318    double Sy, Sx, Sxy, Sx2, Sy2,r_up, r_low, x, y;
1319    double *r;
1320    int a;
1321
1322    r=vcalloc ( 3, sizeof (double));
1323    Sy=Sx=Sxy=Sx2=Sy2=0;
1324
1325    for ( a=0; a<n; a++)
1326      {
1327        x=list[a][0];
1328        y=list[a][1];
1329        Sy+=y;
1330        Sx+=x;
1331        Sy2+=y*y;
1332        Sx2+=x*x;
1333        Sxy+=x*y;
1334      }
1335    r_up=n*Sxy-(Sx*Sy);
1336    r_low=(n*Sx2-(Sx*Sx))*(n*Sy2-(Sy*Sy));
1337    r_low=sqrt(r_low);
1338    r[0]=(r_low==0)?0:r_up/r_low;
1339
1340    x=Sx/(double)n;
1341    y=Sy/(double)n;
1342    r_up=Sxy-(n*x*y);
1343    r_up*=r_up;
1344    r_low=(Sx2-n*x*x)*(Sy2-n*y*y);
1345    r[1]=(r_low==0)?0:r_up/r_low;
1346    r[2]=n;
1347    return r;
1348   }
1349
1350 #define INVERT_LIST(type,wf,rf,function,swap_function)\
1351 type* function (type *list, int len)\
1352     {\
1353     int a, b;\
1354     for ( a=0, b=len-1; a<b; a++, b--)swap_function ( &list[a], &list[b], 1);\
1355     return list;\
1356     }
1357 INVERT_LIST(short,write_size_short,read_size_short, invert_list_short,swap_short)
1358 INVERT_LIST(char,write_size_char,read_size_char,invert_list_char,swap_char)
1359 INVERT_LIST(int,write_size_int,read_size_int,invert_list_int,swap_int)
1360 INVERT_LIST(float,write_size_float,read_size_float,invert_list_float,swap_float)
1361 INVERT_LIST(double,write_size_double,read_size_double,invert_list_double,swap_double)
1362
1363 #define SWAP_FUNCTION(type,wf,rf,function)\
1364 void function(type *a, type *b, int n)\
1365     {\
1366     type t;\
1367     int c;\
1368     for ( c=0;c<n;c++)\
1369         {t=a[c];\
1370          a[c]=b[c];\
1371          b[c]=t;\
1372         }\
1373     }
1374 SWAP_FUNCTION(short,write_size_short,read_size_short,swap_short)
1375 SWAP_FUNCTION(char,write_size_char,read_size_char,swap_char)
1376 SWAP_FUNCTION(int,write_size_int,read_size_int,swap_int)
1377 SWAP_FUNCTION(float,write_size_float,read_size_float,swap_float)
1378 SWAP_FUNCTION(double,write_size_double,read_size_double,swap_double)
1379
1380 #define RETURN_MAX_HORIZ(type,wf,rf,function,comparison,undef)\
1381 type function  (type ** array, int len_array, int field)\
1382     {\
1383     type max;\
1384     int a;\
1385     if ( len_array==0)return 0;\
1386     else\
1387         {\
1388         max=array[field][0];\
1389         for ( a=1; a< len_array; a++)\
1390             if ( array[field][a]!=undef) max=( array[field][a] comparison max)?array[field][a]:max;\
1391         return (int)max;\
1392         }\
1393     }
1394 RETURN_MAX_HORIZ(short,write_size_short,read_size_short,return_max_short_hor,>,UNDEFINED_SHORT)
1395 RETURN_MAX_HORIZ(char,write_size_char,read_size_char,return_max_char_hor,>,UNDEFINED_CHAR)
1396 RETURN_MAX_HORIZ(int,write_size_int,read_size_int,return_max_int_hor,>,UNDEFINED_INT)
1397 RETURN_MAX_HORIZ(float,write_size_float,read_size_float,return_max_float_hor,>,UNDEFINED_FLOAT)
1398 RETURN_MAX_HORIZ(double,write_size_double,read_size_double,return_max_double_hor,>,UNDEFINED_DOUBLE)
1399
1400 RETURN_MAX_HORIZ(short,write_size_short,read_size_short,return_min_short_hor,<,UNDEFINED_SHORT)
1401 RETURN_MAX_HORIZ(char,write_size_char,read_size_char,return_min_char_hor,<,UNDEFINED_CHAR)
1402 RETURN_MAX_HORIZ(int,write_size_int,read_size_int,return_min_int_hor,<,UNDEFINED_INT)
1403 RETURN_MAX_HORIZ(float,write_size_float,read_size_float,return_min_float_hor,<,UNDEFINED_FLOAT)
1404 RETURN_MAX_HORIZ(double,write_size_double,read_size_double,return_min_double_hor,<,UNDEFINED_DOUBLE)
1405
1406
1407
1408 #define BEST_OF_MANY(type,wf,rf,function,undef)\
1409 type function (int n, ...)\
1410         {\
1411         va_list ap;\
1412         int *fop,a;\
1413         type v, best;\
1414         int maximise;\
1415         /*first Arg: number of values\
1416           2nd   Arg: maximise(1)/minimise(0)\
1417           3rd   Arg: *int contains the indice of the best value\
1418           ...   Arg: n type values\
1419         */\
1420         va_start (ap, n);\
1421         maximise=va_arg (ap, int);\
1422         fop=va_arg (ap, int*);\
1423         best=va_arg (ap, type);\
1424         fop[0]=0;\
1425         for ( a=1; a<n; a++)\
1426                 {\
1427                 v=va_arg (ap, type);\
1428                 if (best==undef)\
1429                         {\
1430                         best=v;\
1431                         fop[0]=a;\
1432                         }\
1433                 if ( best==undef || v==undef);\
1434                 else if ( maximise==1 && v>best)\
1435                         {\
1436                         fop[0]=a;\
1437                         best=v;\
1438                         }\
1439                 else if ( maximise==0 && v<best)\
1440                         {\
1441                         fop[0]=a;\
1442                         best=v;\
1443                         }\
1444                 }\
1445         va_end (ap);\
1446         return best;\
1447         }
1448      /*BEST_OF_MANY(short,write_size_short,read_size_short, best_short,UNDEFINED_SHORT)*/
1449      /*BEST_OF_MANY(char,write_size_char,read_size_char,best_char,UNDEFINED_CHAR)*/
1450 BEST_OF_MANY(int,write_size_int,read_size_int,best_int,UNDEFINED_INT)
1451      /*BEST_OF_MANY(float,write_size_float,read_size_float,best_float,UNDEFINED_FLOAT)*/
1452 BEST_OF_MANY(double,write_size_double,read_size_double,best_double,UNDEFINED_DOUBLE)
1453 #define IS_DEFINED(type,function,undef)\
1454 int function(int n, ...)\
1455      {\
1456      int i;\
1457      va_list ap;\
1458 \
1459      va_start(ap, n);\
1460      for ( i=0; i< n; i++)\
1461          {\
1462          if(va_arg(ap,type)==undef)\
1463               {\
1464                 va_end(ap);\
1465                 return 0;\
1466               }\
1467          }\
1468      va_end(ap);\
1469      return 1;\
1470      }
1471      /*IS_DEFINED(short,is_defined_short,UNDEFINED_SHORT)*/
1472      /*IS_DEFINED(char,is_defined_char,  UNDEFINED_CHAR)*/
1473 IS_DEFINED(int,is_defined_int,   UNDEFINED_INT)
1474      /*IS_DEFINED(float,is_defined_float, UNDEFINED_FLOAT)*/
1475 IS_DEFINED(double,is_defined_double,UNDEFINED_DOUBLE)
1476
1477   int return_maxlen ( char ** array, int number)
1478     {
1479     int a;
1480     int max=0;
1481     for ( a=0; a< number; a++)
1482         max=( strlen ( array[a])>max)?strlen ( array[a]):max;
1483
1484     return max;
1485     }
1486
1487
1488 int return_minlen ( char ** array, int number)
1489     {
1490     int a;
1491     int min;
1492
1493     min=strlen( array[0]);
1494     for ( a=1; a< number; a++)
1495         min=( strlen ( array[a])>min)?strlen ( array[a]):min;
1496
1497     return min;
1498     }
1499
1500
1501
1502 float return_mean_diff_float ( float **array, int len, int field,float mean)
1503     {
1504     int a;
1505     float b=0;
1506
1507     for ( a=0; a< len; a++)
1508         {
1509         if ( (mean-array[a][field])!=0)
1510                 b+=sqrt((double)((float) ( mean-array[a][field])*(float)(mean-array[a][field])));
1511         }
1512
1513     return ((float)b/(float)len);
1514     }
1515
1516
1517
1518 void inverse_int ( int**array, int len, int field, int max, int min)
1519     {
1520     int a;
1521     for ( a=0; a< len; a++)
1522         array[a][field]=max-array[a][field]+min;
1523     }
1524 void inverse_float ( float**array, int len, int field, int max, int min)
1525     {
1526     int a;
1527     for ( a=0; a< len; a++)
1528         array[a][field]=max-array[a][field]+min;
1529     }
1530 void inverse_2D_float ( float **array, int start, int len, int start_field, int number_field, float max,float min)
1531     {
1532     int a, b;
1533     for ( a=0; a< start+len; a++)
1534         for ( b=start_field; b< start_field+ number_field; b++)
1535             array[a][b]=max-array[a][b]+min;
1536     }
1537
1538 int max_int (int*i, ...)
1539 {
1540   va_list ap;                                   \
1541   int index, best_value=0, value;
1542   int a=0;
1543  //  expects n values : n, &index, i1, v1, i2, v2...., -1
1544      va_start(ap, i);
1545      while ((index=va_arg(ap,int))!=-1)
1546        {
1547          value=va_arg(ap, int);
1548          if ( a==0 || value>best_value)
1549            {
1550              i[0]=index;
1551              best_value=value;
1552              a=1;
1553            }
1554        }
1555      va_end (ap);
1556      return best_value;
1557 }
1558
1559 /*********************************************************************/
1560 /*                                                                   */
1561 /*                         SHELL INTERFACES                          */
1562 /*                                                                   */
1563 /*                                                                   */
1564 /*********************************************************************/
1565 char* getenv4debug (const char * val)
1566 {
1567   /*efficient mean of getting an environment variable: checks only if one DEBUG is on*/
1568   static int check;
1569
1570   if ( !check)
1571     {
1572
1573       if (getenv       ("DEBUG_BLAST"))check=1;
1574       else if ( getenv ("DEBUG_TREE_COMPARE"))check=1;
1575       else if ( getenv ("DEBUG_MALN"))check=1;
1576       else if ( getenv ("DEBUG_EXTRACT_FROM_PDB"))check=1;
1577       else if ( getenv ("DEBUG_LIBRARY"))check=1;
1578       else if ( getenv ("DEBUG_FUGUE"))check=1;
1579
1580       else if ( getenv ("DEBUG_REFORMAT"))check=1;
1581       else if ( getenv ("DEBUG_RECONCILIATION"))check=1;
1582       else if ( getenv ("DEBUG_TMP_FILE"))check=1;
1583       else if ( getenv ("DEBUG_TREE"))check=1;
1584
1585      else if ( getenv ("DEBUG_SEQ_REFORMAT") && strm (PROGRAM, "SEQ_REFORMAT"))check=2;
1586       else if ( getenv ("DEBUG_TCOFFEE") && strm (PROGRAM, "T-COFFEE"))check=2;
1587       else check=-1;
1588     }
1589
1590   if ( check>0 && strm ( val, "DEBUG_TMP_FILE"))
1591     {
1592       return "1";
1593     }
1594
1595   else if ( check==1)
1596     {
1597       return getenv (val);
1598     }
1599   else if ( check==2)
1600     {
1601       return "1";
1602     }
1603   else
1604     return NULL;
1605 }
1606
1607 char* get_env_variable ( const char *var, int mode)
1608         {
1609             /*mode 0: return NULL if variable not set*/
1610             /*mode 1: crash if variable not set*/
1611             if ( !getenv (var))
1612                {
1613                    if (mode==NO_REPORT)return NULL;
1614                    else if ( mode ==IS_NOT_FATAL)
1615                      {
1616                        fprintf ( stderr, "\nYou must set the variable %s [FATAL]\n", var);
1617                        return NULL;
1618                       }
1619                    else
1620                       {
1621                           fprintf ( stderr, "\nYou must set the variable %s [FATAL]\n", var);
1622                           myexit (EXIT_FAILURE);
1623                           return NULL;
1624                       }
1625                }
1626             else return getenv (var);
1627         }
1628
1629 void get_pwd ( char *name)
1630         {
1631         char *string;
1632         char command[1000];
1633         FILE *fp;
1634
1635
1636         string=vtmpnam(NULL);
1637         sprintf ( command, "pwd > %s", string);
1638         my_system (command);
1639         fp=vfopen ( string, "r");
1640         fscanf ( fp, "%s",name);
1641         vfclose (fp);
1642         sprintf ( command, "rm %s", string);
1643         my_system ( command);
1644         }
1645 int pg_is_installed ( char *pg)
1646         {
1647         char *fname;
1648         char command[1000];
1649         FILE *fp;
1650         int r=0;
1651
1652         return 1;
1653
1654         fname= vtmpnam(NULL);
1655
1656         sprintf ( command, "which %s > %s", pg, fname);
1657         my_system ( command);
1658
1659
1660         if ((fp=find_token_in_file ( fname, NULL, "Command"))){r=1;vfclose(fp);}
1661
1662
1663         return r;
1664
1665         }
1666
1667
1668 /*********************************************************************/
1669 /*                                                                   */
1670 /*                           MISC                                    */
1671 /*                                                                   */
1672 /*********************************************************************/
1673 char *num2plot (int value, int max, int line_len)
1674         {
1675                int   len;
1676                int   value_len;
1677                char *buf;
1678         static char *string;
1679
1680         if ( string==NULL)string=vcalloc (1000, sizeof(char));
1681
1682         if ( line_len==-1)len=30;
1683         else len=line_len;
1684
1685         value_len=((float)value/(float)max)*(float)len;
1686         if ( value==0)
1687             sprintf ( string, "|");
1688         else
1689             {
1690             buf=generate_string(value_len, '*');
1691             sprintf ( string,"%s", buf);
1692             vfree(buf);
1693             }
1694         return string;
1695         }
1696
1697 int   perl_strstr ( char *string, char *pattern)
1698 {
1699   char *tmp;
1700   FILE *fp;
1701   int r;
1702   char command[10000];
1703   char *string2;
1704
1705   if (!string)  return 0;
1706   if (!pattern) return 0;
1707
1708
1709
1710   string2=vcalloc ( strlen (string)+1, sizeof (char));
1711   sprintf ( string2,"%s", string);
1712   string2=substitute (string2, "(", " ");
1713   string2=substitute (string2, ")", " ");
1714   string2=substitute (string2, "'", " ");
1715   tmp=vtmpnam(NULL);
1716   sprintf (command, "perl -e '$s=\"%s\";$x=($s=~/%s/);$x=($x==1)?1:0;print $x;'>%s", string2, pattern,tmp);
1717   system ( command);
1718   if (check_file_exists(tmp))
1719     {
1720       fp=vfopen (tmp, "r");
1721       fscanf (fp, "%d", &r);
1722       vfclose (fp);
1723     }
1724   else
1725     {
1726       fprintf ( stderr, "COM: %s\n");
1727       r=0;
1728     }
1729   vfree (string2);
1730   return r;
1731 }
1732 float grep_function ( char *pattern, char *file)
1733         {
1734         char command [1000];
1735         int a, b, l;
1736         char buf1[100];
1737         char buf2[100];
1738         FILE*fp;
1739         char *s;
1740         float f;
1741
1742
1743
1744         s=vtmpnam(NULL);
1745
1746         sprintf ( command, "grep %s %s > %s",pattern,file, s);
1747         my_system ( command);
1748         if ((fp=vfopen (s, "r"))==NULL )return 0;
1749         else
1750                 {
1751                 fgets ( command, 900, fp);
1752                 l=strlen ( command);
1753                 while ( !isdigit (command[l]))l--;
1754                 a=0;
1755                 while ( isdigit (command[l]) || command[l]=='.')
1756                         {
1757                         buf1[a++]=command[l];
1758                         l--;
1759                         }
1760                 buf1[a]='\0';
1761                 l=strlen (buf1);
1762                 for ( a=0, b=l-1;a< l; a++, b--)
1763                         buf2[b]=buf1[a];
1764                 buf2[l]='\0';
1765
1766                 sscanf ( buf2, "%f", &f);
1767                 sprintf ( command,"rm %s", s);
1768                 my_system ( command);
1769                 vfclose (fp);
1770                 return f;
1771                 }
1772         }
1773
1774 void crash_if ( int val, char *s)
1775     {
1776     if ( val==0)crash(s);
1777     }
1778 void crash ( char *s)
1779         {
1780         int *a;
1781
1782
1783
1784         fprintf ( stderr, "%s",s);
1785         a=vcalloc ( 10, sizeof (int));
1786         a[20]=1;
1787         error_exit();
1788         }
1789
1790 static int *local_table;
1791 int ** make_recursive_combination_table ( int tot_n_param, int *n_param, int *nc, int**table, int field)
1792     {
1793     int a, b, c;
1794
1795     /* makes a table of all possible combinations*/
1796
1797     if ( tot_n_param==0)
1798         {
1799             nc[0]=1;
1800             fprintf ( stderr, "\nNULL RETURNED");
1801             return NULL;
1802         }
1803     if (table==NULL)
1804         {
1805         if ( local_table!=NULL)vfree (local_table);
1806         local_table=vcalloc ( tot_n_param, sizeof (int));
1807         field=0;
1808         for ( a=0; a< tot_n_param; a++)local_table[a]=-1;
1809         for ( a=0; a< tot_n_param; a++)nc[0]=nc[0]*n_param[a];
1810
1811
1812         table=declare_int ( nc[0],tot_n_param);
1813         nc[0]=0;
1814         }
1815
1816     for ( b=0; b<n_param[field]; b++)
1817                {
1818
1819                local_table[field]=b;
1820                if ( field<tot_n_param-1)
1821                   {
1822                   table=make_recursive_combination_table ( tot_n_param, n_param, nc, table, field+1);
1823                   }
1824                else
1825                   {
1826                   for ( c=0; c< tot_n_param; c++)table[nc[0]][c]=local_table[c];
1827                   nc[0]++;
1828                   }
1829                }
1830     return table;
1831     }
1832
1833 /*********************************************************************/
1834 /*                                                                   */
1835 /*                         STRING PROCESSING                         */
1836 /*                                                                   */
1837 /*                                                                   */
1838 /*********************************************************************/
1839 char *strnrchr ( char *s,char x, int n)
1840 {
1841   int a;
1842   for (a=0; a< n; a++)if (s[a]=='x')return s+a;
1843   return NULL;
1844 }
1845 int intlen (int n)
1846 {
1847   char buf [100];
1848   sprintf ( buf, "%d", n);
1849   return strlen (buf)+1;
1850 }
1851 char * update_string (char string1[], char string2[])
1852 {
1853   if ( string1==string2);
1854   else if ( string2==NULL)
1855     {
1856       if ( string1==NULL)string1=vcalloc ( 1, sizeof(char));
1857       string1[0]='\0';
1858     }
1859   else
1860     {
1861       int l1, l2;
1862       l1=read_array_size_new (string1);
1863       l2=strlen (string2)+1;
1864       if (l1<l2)
1865         string1=vrealloc (string1, (l2)*sizeof (char));
1866       sprintf ( string1, "%s", string2);
1867     }
1868   return string1;
1869 }
1870
1871 char* strcatf  (char *string1,char *string2, ...)
1872 {
1873
1874   va_list ap;
1875   char *buf;
1876
1877   buf=vcalloc ( read_array_size_new (string1)+1, sizeof (char));
1878   va_start (ap, string2);
1879   vsprintf (buf, string2, ap);
1880   va_end (ap);
1881   string1=vcat (string1, buf);
1882   vfree (buf);
1883   return string1;
1884 }
1885
1886 char *vcat ( char *st1, char *st2)
1887 {
1888   int l=0;
1889   char *ns;
1890
1891   if ( !st1 && !st2)return NULL;
1892   if ( st1) l+=strlen (st1);
1893   if ( st2) l+=strlen (st2);
1894   l++;
1895
1896   ns=vcalloc ( l, sizeof (char));
1897   sprintf ( ns, "%s%s", (st1)?st1:"", (st2)?st2:"");
1898   return ns;
1899 }
1900 int print_param ( char *param, FILE *fp);
1901 int print_param ( char *param, FILE *fp)
1902 {
1903   static char **p;
1904   static int np;
1905   int a;
1906
1907   if (!p)p=vcalloc (1000, sizeof (char*));
1908
1909   for ( a=0; a<np; a++)if (p[a]==param) return 0;
1910
1911   p[np++]=param;
1912
1913   fprintf ( fp, "\nPARAMETERS: %s\n", param);
1914   return 1;
1915 }
1916
1917 int strget_param ( char *string, char *token1, char *token2, char *format, ...)
1918 {
1919   /*string: command line
1920     token1: parameter
1921     token2: default value (in a string)
1922     format: standard scanf format
1923     ... arguments for scanf
1924   */
1925   char *buf;
1926   char *buf2;
1927   int n;
1928   va_list ap;
1929
1930   if (!string) return 0;
1931   //print_param (string, stdout);
1932   //print_param (string, stderr);
1933   va_start (ap, format);
1934   buf=after_strstr ( string, token1);
1935
1936
1937   if (buf)
1938     {
1939       buf2=vcalloc (strlen (buf)+1, sizeof (char));
1940       sprintf ( buf2, "%s", buf);
1941       buf2=substitute (buf2, "__", " ");
1942       n=my_vsscanf (buf2, format, ap);
1943       vfree (buf2);
1944     }
1945   else     {n=my_vsscanf (token2, format, ap);}
1946   va_end (ap);
1947
1948   return n;
1949 }
1950
1951 int strscanf (char *string1,char *token, char *format,...)
1952 {
1953   char *buf;
1954   va_list ap;
1955   int n;
1956
1957   va_start (ap,format);
1958   buf=after_strstr ( string1, token);
1959   if ( buf){n=my_vsscanf (buf, format, ap);}
1960   else n=0;
1961
1962   va_end (ap);
1963   return n;
1964 }
1965
1966 int match_motif ( char *string, char **motif)//crude string matching, the motif and the string have the same length
1967   {
1968   int l, a;
1969   if (!motif) return 0;
1970   l=strlen (string);
1971   for ( a=0; a<l; a++)
1972       {
1973         if ( motif[a][0]!='*' && !strchr (motif[a], string[a]))return 0;
1974       }
1975   return 1;
1976 }
1977
1978 char *after_strstr ( char *string, char *token)
1979 {
1980   char *p;
1981   if ( (p=vstrstr (string, token)))return p+strlen (token);
1982   else return NULL;
1983 }
1984 char* lstrstr ( char *in, char *token)//return strstr if matches on the left
1985 {
1986   char *s;
1987   if ((s=vstrstr(in, token))&&s==in)return s;
1988   else return NULL;
1989
1990 }
1991
1992 char *vstrstr ( char *in, char *token)
1993 {
1994   if (!in || !token) return NULL;
1995   else return strstr (in, token);
1996 }
1997 char ** push_string (char *val, char **stack, int *nval, int position)
1998 {
1999   char **new_stack;
2000   int a;
2001
2002
2003   if (!val || nval[0]<=position)return stack;
2004   nval[0]++;
2005   new_stack=vcalloc ( nval[0], sizeof (char*));
2006   new_stack[position]=val;
2007   for (a=0; a< position; a++)new_stack[a]=stack[a];
2008   for (a=position+1; a<nval[0]; a++)new_stack[a]=stack[a-1];
2009   vfree (stack);
2010
2011   return new_stack;
2012 }
2013
2014 int vsrand (int val)
2015 {
2016   static int initialized;
2017
2018   if (initialized) return 0;
2019   else if (   getenv ("DEBUG_SRAND"))
2020     {
2021       srand (10);
2022       initialized=1;
2023     }
2024   else
2025     {
2026
2027       /*Make sure two processes in a row do not get the same time value*/
2028       srand ((val==0)?time(NULL):val);
2029       initialized=1;
2030     }
2031   return 1;
2032 }
2033 int  *randomize_list (int *list, int len, int ncycle)
2034 {
2035   int p1, p2, a, buf;
2036
2037   vsrand (0);
2038   if ( ncycle==0)ncycle=len;
2039   for ( a=0; a<ncycle; a++)
2040     {
2041       p1=rand()%len;
2042       p2=rand()%len;
2043       buf=list[p1];
2044       list[p1]=list[p2];
2045       list[p2]=buf;
2046     }
2047   return list;
2048 }
2049 /*Replace by a gap the parts of the two strings that do not OVERLAP*/
2050 /*returns the length of the overlap*/
2051 int vstrcmp (const char *s1, const char *s2)
2052 {
2053   if ( !s1 && !s2)return 0;
2054   else if ( !s1 || !s2)return 1;
2055   else return strcmp ( s1, s2);
2056 }
2057 int vstrncmp (const char *s1, const char *s2, int n)
2058 {
2059   if ( !s1 && !s2)return 0;
2060   else if ( !s1 || !s2)return 1;
2061   else return strncmp ( s1, s2, n);
2062 }
2063 FILE *print_array_char (FILE *out, char **array, int n, char *sep)
2064         {
2065         int a;
2066         if ( array==NULL || read_size_char (array,sizeof (char*))<n)
2067            {
2068            fprintf ( stderr, "\nORB in print_array_char [FATAL]\n");
2069            crash("");
2070            }
2071         for ( a=0; a< n; a++)fprintf ( out, "%s%s", array[a],sep);
2072         return out;
2073         }
2074 char * path2filename ( char *array)
2075 {
2076   Fname *F;
2077   char *name;
2078
2079   name=vcalloc ( strlen (array)+2, sizeof (char));
2080   F=parse_fname (array);
2081   if ( F->suffix)sprintf ( name, "%s.%s", F->name, F->suffix);
2082   else sprintf (name, "%s", F->name);
2083   free_fname (F);
2084   return name;
2085 }
2086 Fname* parse_fname ( char *array)
2087          {
2088          int l;
2089          Fname *F;
2090
2091
2092
2093          F=declare_fname (sizeof (array));
2094
2095          sprintf ( F->full, "%s", array);
2096          sprintf ( F->path, "%s", array);
2097          l=strlen (array);
2098          while (l!=-1 && (F->path)[l]!='/')(F->path)[l--]='\0';
2099
2100          sprintf ( F->name, "%s", array+l+1);
2101          l=strlen (F->name);
2102          while (l!=-1)
2103            {
2104             if((F->name)[l]=='.')
2105               {
2106               F->name[l]='\0';
2107               sprintf ( F->suffix, "%s", F->name+l+1);
2108               break;
2109               }
2110             else l--;
2111             }
2112
2113         return F;
2114         }
2115 char *filename2path (char *name)
2116 {
2117   char *nname;
2118   int x;
2119   if (isdir (name))return name;
2120
2121   x=strlen (name)-1;
2122   nname=vcalloc (x+2, sizeof (char));
2123   sprintf ( nname, "%s", name);
2124   while ( x >=0 && nname[x]!='/')nname[x--]='\0';
2125
2126   if ( !isdir (nname) || !nname[0]){vfree (nname); return NULL;}
2127   return nname;
2128 }
2129
2130
2131
2132
2133
2134 char *extract_suffixe ( char *array)
2135         {
2136         int l;
2137         char *new_string;
2138         char *x;
2139         l=strlen (array);
2140         new_string=vcalloc ( l+1, sizeof (char));
2141         sprintf (new_string, "%s",array);
2142
2143         x=new_string+l;
2144         while (x!=new_string && x[0]!='.' && x[0]!='/' )x--;
2145         if ( x[0]=='.')x[0]='\0';
2146         else if (x[0]=='/')return x+1;
2147
2148         while ( x!=new_string && x[0]!='/')x--;
2149
2150         return (x[0]=='/')?x+1:x;
2151         }
2152 void string_array_upper ( char **string, int n)
2153      {
2154      int a;
2155      for ( a=0; a< n; a++)upper_string (string[a]);
2156      }
2157 void string_array_lower ( char **string, int n)
2158      {
2159      int a;
2160      for ( a=0; a< n; a++)lower_string (string[a]);
2161      }
2162
2163 char *upper_string ( char *string)
2164         {
2165         int len, a;
2166
2167         len=strlen ( string);
2168         for ( a=0; a< len; a++)string[a]=toupper ( string[a]);
2169         return string;
2170         }
2171 char *lower_string ( char *string)
2172         {
2173         int len, a;
2174
2175         len=strlen ( string);
2176         for ( a=0; a< len; a++)string[a]=tolower ( string[a]);
2177         return string;
2178         }
2179 void string_array_convert ( char **array, int n_strings, int ns, char **sl)
2180         {
2181         int a;
2182
2183         for ( a=0; a< n_strings; a++)string_convert ( array[a], ns, sl);
2184         }
2185 void string_convert( char *string, int ns, char **sl)
2186         {
2187         int a, l;
2188         l=strlen ( string);
2189         for ( a=0; a< l; a++)
2190             string[a]=convert(string[a], ns, sl);
2191         }
2192 int convert ( char c, int ns, char **sl)
2193         {
2194         int a;
2195         int return_char;
2196
2197         for ( a=0; a< ns; a++)
2198             {
2199             if ((return_char=convert2 ( c, sl[a]))!=-1)
2200                 return return_char;
2201             }
2202         return c;
2203
2204
2205         }
2206 int convert2 ( char c, char *list)
2207     {
2208     int a;
2209     int l1;
2210     int return_char;
2211
2212     l1=strlen ( list);
2213
2214     return_char=(list[l1-1]=='#')?c:list[l1-1];
2215
2216     for ( a=0; a< l1; a++)
2217            if (list[a]=='#')return return_char;
2218            else if ( list[a]==c)return return_char;
2219
2220     return -1;
2221     }
2222 char* substitute_old ( char *string_in, char *t, char *r)
2223 {
2224   char *string_out;
2225   char *p, *heap_in;
2226   int delta, l;
2227   /*REplaces every occurence of token t with token r in string_in*/
2228
2229   if ( string_in==NULL || t==NULL || r==NULL) return string_in;
2230
2231   heap_in=string_in;
2232
2233   l=read_array_size_new ((void*)string_in)+1;
2234
2235   string_out=vcalloc (l, sizeof (char));
2236   delta=strlen(r)-strlen (t);
2237   delta=(delta<0)?0:delta;
2238
2239   while ( (p=strstr ( string_in, t))!=NULL)
2240     {
2241
2242       p[0]='\0';
2243       if ( delta)
2244         {
2245           l+=delta;
2246           string_out=vrealloc(string_out, sizeof (char)*l);
2247         }
2248
2249       strcat ( string_out, string_in);
2250       strcat ( string_out, r);
2251       string_in=p+strlen (t);
2252     }
2253   strcat ( string_out, string_in);
2254   if (l<strlen (string_out))
2255     {
2256       heap_in=vrealloc (heap_in, sizeof(char)*(strlen (string_out)+1));
2257     }
2258   sprintf ( heap_in, "%s", string_out);
2259   vfree (string_out);
2260   return heap_in;
2261 }
2262 //Makes sure substitutions of the tild only occur on valid UNIX pathnames
2263 char* tild_substitute ( char *string_in, char *t, char *r)
2264 {
2265   char *p;
2266
2267   p=strstr ( string_in, "~");
2268   if ( p==NULL)return string_in;
2269   else if (p && p!=string_in && p[-1]!='/')return string_in;
2270
2271   else
2272     {
2273       return substituteN (string_in, t, r,1);
2274     }
2275
2276 }
2277 char* substitute_char ( char *string_in, char t, char r)
2278 {
2279   int n=0;
2280   while ( string_in && string_in[n]!='\0')
2281     {
2282       if ( string_in[n] == t)string_in[n]=r;
2283       n++;
2284     }
2285   return string_in;
2286 }
2287 char* substitute ( char *string_in, char *t, char *r)
2288 {
2289   /*REplaces every occurence of token t with token r in string_in*/
2290   return substituteN(string_in, t, r,0);
2291 }
2292 char* substituteN ( char *string_in, char *t, char *r, int N)
2293 {
2294   char *string_out;
2295   char *p, *heap_in, n=0;
2296   int delta, l, lsi,lso,lr,lt,nt;
2297
2298   /*REplaces the first N Occurences*/
2299   if ( string_in==NULL || t==NULL || r==NULL) return string_in;
2300
2301   heap_in=string_in;
2302
2303   l=read_array_size_new ((void*)string_in)+1;
2304   lr=strlen(r);
2305   lt=strlen(t);
2306   lsi=strlen (string_in);
2307
2308   delta=((lr-lt)>0)?(lr-lt):0;
2309   nt=0;
2310   while ( (p=strstr (string_in, t))!=NULL)
2311     {
2312       string_in=p+lt;
2313       nt++;
2314     }
2315   string_in=heap_in;
2316
2317   lso=nt*delta+lsi;
2318   string_out=vcalloc (lso+1, sizeof (char));
2319
2320   while ((N==0 ||n<N) && (p=strstr ( string_in, t))!=NULL)
2321     {
2322       p[0]='\0';
2323       strcat ( string_out, string_in);
2324       strcat ( string_out, r);
2325       string_in=p+lt;
2326       if (N!=0)n++;
2327     }
2328   strcat ( string_out, string_in);
2329   if (l<(lso+1)){HERE ("Realloc %s %s delta=%d %d %d", t, r,delta, lsi, lso);heap_in=vrealloc (heap_in, sizeof (char)*(lso +1));}
2330   sprintf ( heap_in, "%s", string_out);
2331   vfree (string_out);
2332   return heap_in;
2333 }
2334
2335
2336 char **clean_string (int n, char **string)
2337 {
2338   int a,b,c,l;
2339   if ( !string) return string;
2340   for (a=0; a< n; a++)
2341     {
2342       if (!string[a]) continue;
2343       l=strlen (string[a]);
2344       for (b=0; b<l; b++)
2345         {
2346           c=string[a][b];
2347           if (!isgraph(c) && c!='\n' && c!='\t' && c!=' ')string[a][b]=' ';
2348         }
2349     }
2350   return string;
2351 }
2352
2353
2354
2355
2356 int str_overlap ( char *string1, char *string2, char x)
2357         {
2358         int a, b;
2359         int l1, l2;
2360         int **array;
2361         int max1=0, max2=0;
2362         int score;
2363         int max_val=-1;
2364         int end1, end2;
2365         int start1, start2;
2366
2367         if ( strm ( string1, string2))return 0;
2368         else
2369             {
2370             l1=strlen ( string1);
2371             l2=strlen ( string2);
2372             array=declare_int ( strlen ( string1), strlen ( string2));
2373             for ( a=0; a< l1; a++)
2374                 for ( b=0; b< l2; b++)
2375                     {
2376                         if ( a==0 || b==0)array[a][b]=(string1[a]==string2[b])?1:0;
2377                         else
2378                             if (string1[a]==string2[b])
2379                                {
2380                                score=array[a][b]=array[a-1][b-1]+1;
2381                                if ( max_val<score)
2382                                   {
2383                                   max_val=score;
2384                                   max1=a;
2385                                   max2=b;
2386                                   }
2387                                }
2388                     }
2389             start1=(max1+1)-max_val;
2390             end1=  max1;
2391
2392             start2=(max2+1)-max_val;
2393             end2=  max2;
2394
2395             for ( a=0; a< l1; a++)if ( a<start1 || a> end1)string1[a]=x;
2396             for ( a=0; a< l2; a++)if ( a<start2 || a> end2)string2[a]=x;
2397
2398             free_int ( array, l1);
2399
2400             return max_val;
2401             }
2402         }
2403
2404 int get_string_line ( int start, int n_lines, char *in, char *out)
2405         {
2406         int nl=0;
2407         int a=0;
2408         int c=0;
2409
2410         while ( nl<n_lines)
2411                 {
2412                 while ( (c=in[start++])!='\n' && c!='\0')
2413                         {
2414                         out[a++]=c;
2415                         }
2416                 out[a++]='\n';
2417                 nl++;
2418                 }
2419         out[a]='\0';
2420         return (c=='\0')?-1:start;
2421         }
2422
2423
2424 FILE * output_string_wrap ( int wrap,char *string, FILE *fp)
2425         {
2426
2427          int a, b,l;
2428
2429          l=strlen ( string);
2430
2431          for ( a=0, b=1; a< l; a++, b++)
2432                 {
2433                 fprintf ( fp, "%c", string[a]);
2434                 if ( b==wrap)
2435                         {
2436                         fprintf ( fp, "\n");
2437                         b=0;
2438                         }
2439                 }
2440          return fp;
2441          }
2442
2443 char * extract_char ( char * array, int first, int len)
2444     {
2445     char *array2;
2446     int a;
2447
2448     len= ( len<0)?0:len;
2449
2450     array2=vcalloc ((len+1), sizeof (char));
2451
2452     for ( a=0; a<len; a++)
2453         array2[a]=array[first++];
2454
2455     array2[a]='\0';
2456
2457     return array2;
2458     }
2459
2460 int check_cl4t_coffee (int argc, char **argv)
2461 {
2462   char *name;
2463   if ( name_is_in_list ("-other_pg", argv, argc, 100)!=-1)return 1;
2464   else if (name_is_in_list ( "tcoffee_test_seq.pep", argv, argc, 100)!=-1)return 1;
2465   else
2466     {
2467       int a,  inseq, result;
2468       FILE *fp;
2469       char command[10000];
2470       command[0]='\0';
2471
2472       for ( inseq=0,a=0; a<argc; a++)
2473         {
2474           if ( inseq && argv[a][0]=='-')
2475             {
2476               inseq=0;
2477             }
2478           else if (strm ( argv[a], "-seq"))inseq=1;
2479
2480           if ( inseq==0)
2481             {
2482               strcat ( command, " ");
2483               strcat ( command, argv[a]);
2484             }
2485         }
2486       name=vcalloc ( 100, sizeof (char));
2487       sprintf ( name, "TCIT_%d", (int)rand()%10000);
2488       strcat ( command, " -seq tcoffee_test_seq.pep -quiet -no_error_report");
2489       fp=vfopen ( name, "w");
2490       fprintf ( fp, ">A\nthecat\n>B\nthecat\n");
2491       vfclose (fp);
2492       result=safe_system (command);
2493       printf_system ( "rm %s.*", name);
2494       vfree (name);
2495       if (result) {myexit (EXIT_FAILURE);return 0;}
2496       else return 1;
2497     }
2498 }
2499
2500 char** merge_list ( char **argv, int *argc)
2501        {
2502        int a, b;
2503        int n_in;
2504        char **out;
2505        char current [STRING];
2506
2507        out=declare_char (argc[0], STRING);
2508        n_in=argc[0];
2509        argc[0]=0;
2510
2511        a=0;
2512        while (a< n_in && !is_parameter ( argv[a]))
2513          {
2514            sprintf (out[argc[0]++], "%s",  argv[a]);
2515            argv[a][0]='\0';
2516            a++;
2517          }
2518
2519
2520        for ( a=0; a< n_in; a++)
2521          {
2522            if ( is_parameter (argv[a]))
2523              {
2524                sprintf ( out[argc[0]++], "%s", argv[a]);
2525                sprintf ( current, "%s", argv[a]);
2526
2527                for ( b=0; b< n_in;)
2528                  {
2529                    if ( is_parameter (argv[b]) && strm (current, argv[b]))
2530                         {
2531                           argv[b][0]='\0';
2532                           b++;
2533                           while ( b<n_in && !is_parameter ( argv[b]) )
2534                             {
2535                               if (argv[b][0])
2536                                {
2537                                  sprintf ( out[argc[0]++], "%s", argv[b]);
2538                                  argv[b][0]='\0';
2539                                }
2540                              b++;
2541                             }
2542                         }
2543
2544                    else b++;
2545
2546
2547
2548                  }
2549              }
2550          }
2551
2552        free_char (argv, -1);
2553        return out;
2554        }
2555
2556
2557 int *  string2num_list_old ( char *string)
2558 {
2559   /*Breaks down a list of numbers separated by any legal separator and put them in an array of the right size*/
2560   /*Returns list a list of integer*/
2561   /*Skips non numbers*/
2562
2563   char *buf, *s;
2564   int  *list, n;
2565
2566
2567   buf=vcalloc ( strlen (string)+1, sizeof (char));
2568
2569   n=0;
2570   sprintf ( buf, "%s", string);
2571   s=strtok (buf, SEPARATORS);
2572   while (s!=NULL)
2573           {
2574             n++;
2575             s=strtok (NULL, SEPARATORS);
2576           }
2577   list=vcalloc (n+1, sizeof (int));
2578
2579   n=0;
2580   sprintf ( buf, "%s", string);
2581   s=strtok (buf, SEPARATORS);
2582   while (s!=NULL)
2583           {
2584             if (is_number(s))list[n++]=atoi(s);
2585             s=strtok (NULL, SEPARATORS);
2586           }
2587   vfree (buf);
2588
2589   return list;
2590 }
2591
2592
2593
2594 int *name_array2index_array ( char **list1, int n1, char **list2, int n2)
2595 {
2596   int *list,a, max;
2597   /*returns an indexed list of the position of list1 in list2*/
2598   list=vcalloc ( n1, sizeof (int));
2599   for ( a=0, max=0; a< n1; a++)max=MAX(max, strlen (list1[a]));
2600   for ( a=0       ; a< n2; a++)max=MAX(max, strlen (list2[a]));
2601
2602   for ( a=0; a< n1; a++)
2603     {
2604       list[a]=name_is_in_list (list1[a],list2,n2,max);
2605     }
2606   return list;
2607 }
2608
2609 int * string2num_list ( char *string)
2610 {
2611   return string2num_list2(string, SEPARATORS);
2612 }
2613 int * string2num_list2 ( char *string, char *separators)
2614 {
2615    /*Breaks down a list of numbers separated by any legal separator and put them in an array of the right size*/
2616   /*Returns list a list of integer*/
2617   /*Skips non numbers*/
2618
2619   char **nlist;
2620   int *list;
2621   int a, m, n;
2622
2623   nlist=string2list2 (string, separators);
2624
2625   if (nlist==NULL) return NULL;
2626
2627   n=atoi(nlist[0]);
2628   list=vcalloc ( n, sizeof (int));
2629
2630   for (m=1, a=1; a< n; a++)
2631     {
2632       if (is_number(nlist[a]))list[m++]=atoi(nlist[a]);
2633     }
2634
2635   list[0]=m;
2636   free_char (nlist, -1);
2637   if ( m==0){vfree(list); return NULL;}
2638   else
2639     {
2640       return list;
2641     }
2642
2643   return NULL;
2644
2645 }
2646 char * list2string  ( char **list, int n)
2647 {
2648   return list2string2 ( list, n, " ");
2649 }
2650 char * list2string2 ( char **list,int n, char* sep)
2651 {
2652   int l, a;
2653   char *string;
2654
2655   for ( l=0,a=0; a<n; a++)l+=strlen ( list[a])+1;
2656   string=vcalloc (l+1, sizeof (char));
2657   for ( a=0; a< n; a++)
2658       {
2659         strcat ( string, list[a]);
2660         strcat ( string, sep);
2661       }
2662     return string;
2663   }
2664
2665 char **  string2list ( char *string)
2666   {
2667     return string2list2(string, SEPARATORS);
2668   }
2669 char **  string2list2 ( char *string, char *separators)
2670   {
2671   /*Breaks down a list of words separated by any legal separator and put them in an array of the right size*/
2672   /*Returns list a list of char with the size written in list[0]*/
2673
2674
2675   char *buf, *s;
2676   char  **list;
2677   int n, max_len;
2678
2679   if ( string==NULL)return NULL;
2680   buf=vcalloc ( strlen (string)+1, sizeof (char));
2681
2682
2683   n=max_len=0;
2684   sprintf ( buf, "%s", string);
2685   s=strtok (buf, separators);
2686
2687   while (s!=NULL)
2688           {
2689             n++;
2690             max_len=MAX(max_len,(strlen (s)));
2691             s=strtok (NULL, separators);
2692
2693           }
2694
2695   if ( n==0){vfree(buf); return NULL;}
2696
2697   list=declare_arrayN (2, sizeof (char), n+1, max_len+1);
2698
2699   n=1;
2700   sprintf ( buf, "%s", string);
2701   s=strtok (buf, separators);
2702   while (s!=NULL)
2703           {
2704             sprintf (list[n++], "%s",s);
2705             s=strtok (NULL, separators);
2706           }
2707
2708   sprintf (list[0], "%d", n);
2709
2710   vfree (buf);
2711
2712   return list;
2713 }
2714
2715 char** break_list ( char **argv, int *argc, char *separators)
2716        {
2717        int a, b;
2718        int n_in;
2719        char **out;
2720        char **ar=NULL;
2721        int n_ar;
2722        int cont=1;
2723
2724        /*Breaks down the argv command line in smaller units, breaking at every separator*/
2725        out=vcalloc (MAX_N_PARAM, sizeof (char*));
2726        n_in=argc[0];
2727        argc[0]=0;
2728
2729        if ( n_in>=MAX_N_PARAM)
2730          {
2731            fprintf ( stderr, "\nERROR: too many parameters, recompile with MAX_N_PARAM set at a higher velue [FATAL:%s]\n", PROGRAM);\
2732            myexit (EXIT_FAILURE);
2733          }
2734
2735        for ( a=0; a< n_in; a++)
2736            {
2737
2738
2739
2740              if (cont)ar=get_list_of_tokens( argv[a], separators,&n_ar);
2741              else ar=get_list_of_tokens( argv[a],"",&n_ar);
2742
2743
2744              for ( b=0; b< n_ar; b++)
2745                {
2746                  out[argc[0]]=vcalloc( strlen (ar[b])+1, sizeof (char));
2747                  sprintf (out[argc[0]++], "%s", ar[b]);
2748                }
2749              free_char (ar, -1);
2750              ar=NULL;
2751              if ( strstr (argv[a], "-other_pg"))cont=0;
2752            }
2753        free_char (ar, -1);
2754        return out;
2755        }
2756
2757 char *invert_string2 (char *string)
2758 {
2759   char *buf;
2760   int a, b, l;
2761
2762   l=strlen (string);
2763   buf=vcalloc ( l+1, sizeof (char));
2764   for ( a=l-1, b=0; a>=0; a--, b++)
2765     buf[b]=string[a];
2766   sprintf (string, "%s", buf);
2767   vfree (buf);
2768   return string;
2769 }
2770 char *invert_string (char *string)
2771 {
2772   return string2inverted_string(string);
2773 }
2774 char* string2inverted_string(char *string)
2775 {
2776   char *buf;
2777   int a, b, l;
2778
2779   l=strlen (string);
2780   buf=vcalloc ( l+1, sizeof (char));
2781   for ( a=l-1, b=0; a>=0; a--, b++)
2782     buf[b]=string[a];
2783   return buf;
2784 }
2785
2786 char ** get_list_of_tokens ( char *in_string, char *separators, int *n_tokens)
2787 {
2788     char **list=NULL;
2789     char *p=NULL;
2790     char *string;
2791
2792
2793     n_tokens[0]=0;
2794     if ( in_string==NULL || strm(in_string, ""));
2795     else if ( in_string[0]=='[')
2796       {
2797         list=declare_char (1, strlen ( in_string)+1);
2798         sprintf ( list[n_tokens[0]], "%s",in_string);
2799         n_tokens[0]++;
2800       }
2801     else
2802       {
2803         list=declare_char (strlen ( in_string)+1, 1);
2804         string=vcalloc ( strlen(in_string)+1, sizeof (char));
2805         sprintf ( string, "%s", in_string);
2806
2807         while ( (p=strtok ((p==NULL)?string:NULL, ((separators==NULL)?SEPARATORS:separators)))!=NULL)
2808            {
2809            list[n_tokens[0]]=vrealloc ( list[n_tokens[0]], sizeof (char) *strlen (p)+1);
2810            sprintf ( list[n_tokens[0]], "%s", p);
2811            n_tokens[0]++;
2812            }
2813
2814         vfree (string);
2815         }
2816    return list;
2817    }
2818
2819 char **ungap_array ( char **array, int n)
2820         {
2821         int a;
2822         for ( a=0; a< n; a++)ungap(array[a]);
2823         return array;
2824         }
2825
2826 void ungap ( char *seq)
2827 {
2828   remove_charset ( seq, "ungap");
2829 }
2830 int seq2len (char *seq, char *pset,char *nset)
2831 {
2832   int a, l, t=0;
2833   //count all the residues in pset and NOT in nset
2834   if ( !seq) return 0;
2835
2836   l=strlen (seq);
2837   //returns the len of the string
2838   for (a=0; a< l; a++)
2839     {
2840       char c=seq[a];
2841       if ( pset && nset && strchr (pset, c) && !strchr (nset, c))t++;
2842       else if ( pset && strchr  (pset, c))t++;
2843       else if ( nset && !strchr (nset, c))t++;
2844     }
2845   return t;
2846 }
2847 int seq2res_len (char *seq)
2848 {
2849   return seq2len (seq, NULL, GAP_LIST);
2850 }
2851 char* remove_charset_from_file (char *fname, char *set)
2852 {
2853   char *tmp;
2854   char c;
2855   FILE *fp1;
2856   FILE *fp2;
2857
2858   fp1=vfopen (fname, "r");
2859   fp2=vfopen (tmp=vtmpnam (NULL), "w");
2860   while ( (c=fgetc(fp1))!=EOF)
2861     {
2862     if (!strchr ( set,c))fprintf ( fp2, "%c", c);
2863     }
2864   vfclose (fp1);
2865   vfclose (fp2);
2866   return tmp;
2867 }
2868
2869 void remove_charset ( char *seq, char *set)
2870         {
2871         int a, b, l;
2872         char *set2;
2873
2874         set2=vcalloc (256, sizeof (char));
2875         if ( strm (set, "!alnum"))
2876           {
2877             for ( b=0,a=1;a< 256; a++)if ( !isalnum (a))set2[b++]=a;
2878           }
2879         else if ( strm ( set, "ungap"))
2880           {
2881             sprintf ( set2, "%s", GAP_LIST);
2882           }
2883         else
2884           {
2885             sprintf ( set2, "%s", set);
2886           }
2887
2888         l=strlen ( seq);
2889         for (b=0, a=0; a<=l; a++)
2890           {
2891             if ( strchr ( set2, seq[a]));
2892             else seq[b++]=seq[a];
2893           }
2894         seq[b]='\0';
2895         vfree (set2);
2896         }
2897
2898
2899 char **char_array2number ( char ** array, int n)
2900        {
2901        int a;
2902        for ( a=0; a< n; a++)array[a]=char2number(array[a]);
2903        return array;
2904        }
2905 char *char2number ( char * array)
2906        {
2907        int a, l;
2908
2909
2910        l=strlen ( array);
2911        for ( a=0; a< l; a++)
2912            {
2913              if ( isdigit(array[a]) && array[a]!=NO_COLOR_RESIDUE && array[a]!=NO_COLOR_GAP )array[a]-='0';
2914              else if ( array[a]<9);
2915              else if ( array[a]==NO_COLOR_RESIDUE || array[a]==NO_COLOR_GAP)array[a]=NO_COLOR_RESIDUE;
2916            }
2917        return array;
2918        }
2919 long atop (char*p)
2920 {
2921   /*turns a char into a pointer*/
2922   if ( p==NULL) return 0;
2923   else return atol(p);
2924 }
2925
2926 char *mark_internal_gaps(char *seq, char symbol)
2927 {
2928   int l, a, gap;
2929   int in_seq;
2930   char *cache_seq;
2931
2932   l=strlen(seq);
2933   cache_seq=vcalloc ( l+1, sizeof (char));
2934   sprintf ( cache_seq, "%s", seq);
2935
2936   for ( gap=0, in_seq=0,a=0; a< l; a++)
2937     {
2938       gap=is_gap(seq[a]);
2939       if ( !gap && !in_seq)in_seq=1;
2940       if (gap && in_seq)seq[a]=symbol;
2941     }
2942
2943   for (gap=0, in_seq=0,a=l-1; a>=0; a--)
2944     {
2945       gap=is_gap(seq[a]);
2946       if ( !gap && !in_seq)break;
2947       if (gap && !in_seq)seq[a]=cache_seq[a];
2948     }
2949   vfree(cache_seq);
2950   return seq;
2951 }
2952
2953 void splice_out ( char *seq, char x)
2954
2955         {
2956         int a, b, l;
2957
2958         l=strlen ( seq);
2959         for (b=0, a=0; a<=l; a++)
2960                 if ( seq[a]==x);
2961                 else seq[b++]=seq[a];
2962         seq[b]='\0';
2963         }
2964 char *splice_out_seg ( char *seq, int pos, int len)
2965 {
2966   int l, a;
2967
2968   if (seq==NULL || pos<0) return seq;
2969   l=strlen (seq);
2970   if ( l<(pos+len))
2971     printf_exit ( EXIT_FAILURE, stderr, "Splice_out_seg out of bound: Length %d seg: [%d %d] [splice_out_seg::util.c][FATAL:%s]\n", l, pos, pos+len, PROGRAM);
2972   l-=len;
2973   for (a=pos; a< l; a++)
2974     seq[a]=seq[a+len];
2975   seq[a]='\0';
2976   return seq;
2977 }
2978
2979 int isblanc ( char *buf)
2980    {
2981    int a, l;
2982
2983    if ( buf==NULL)return 0;
2984    l=strlen (buf);
2985    for ( a=0; a< l; a++)
2986         if (isalnum (buf[a]))return 0;
2987    return 1;
2988    }
2989
2990
2991
2992 int is_number ( char *num)
2993    {
2994      int a, l;
2995      l=strlen (num);
2996
2997      for (a=0;a<l; a++)
2998        if ( !strchr ("0123456789.-+", num[a]))return 0;
2999      return 1;
3000    }
3001
3002 int is_alnum_line ( char *buf)
3003    {
3004    int a, l;
3005    l=strlen (buf);
3006    for ( a=0; a< l; a++)
3007         if (isalnum (buf[a]))return 1;
3008    return 0;
3009    }
3010 int is_alpha_line ( char *buf)
3011    {
3012    int a, l;
3013    l=strlen (buf);
3014    for ( a=0; a< l; a++)
3015         if (isalpha (buf[a]))return 1;
3016    return 0;
3017    }
3018 int case_insensitive_strcmp ( char *string1, char *string2)
3019 {
3020   int a;
3021   int l;
3022
3023   if ( !string1 && !string2) return 1;
3024   else if ( !string1 && !string2) return 0;
3025   else
3026     {
3027       l=strlen (string1);
3028       for ( a=0; a< l; a++)
3029         {
3030           if (tolower(string1[a])!=tolower(string2[a]))return 0;
3031         }
3032     }
3033   return 1;
3034 }
3035 int get_string_sim ( char *string1, char *string2, char *ignore)
3036         {
3037         int len1;
3038         int len2;
3039         int a;
3040         int pos=0;
3041         int sim=0;
3042         char r1, r2;
3043
3044
3045         len1=strlen (string1);
3046         len2=strlen (string2);
3047
3048         if ( len1!=len2)return 0;
3049
3050         for ( a=0; a< len1; a++)
3051                 {
3052                 r1=string1[a];
3053                 r2=string2[a];
3054                 if ( !is_in_set (r1, ignore) && !is_in_set (r2, ignore))
3055                         {
3056                         pos++;
3057                         if (r1==r2)
3058                                 {
3059                                 sim++;
3060                                 }
3061                         }
3062                 }
3063
3064         if (pos==0)
3065                  return 0;
3066         else
3067                 return (int) (sim*100)/pos;
3068
3069         }
3070 int is_aa  ( char x)
3071          {
3072          return (is_in_set (x, AA_ALPHABET) && !is_in_set (x, GAP_LIST));
3073          }
3074 int is_rna ( char x)
3075          {
3076          return (is_in_set (x, RNAONLY_ALPHABET)&& !is_in_set (x, GAP_LIST));
3077          }
3078 int is_dna ( char x)
3079          {
3080          return (is_in_set (x, DNA_ALPHABET)&& !is_in_set (x, GAP_LIST));
3081          }
3082 int is_gap ( char x)
3083          {
3084          return ( is_in_set ( x, GAP_LIST));
3085          }
3086 int is_gop ( int p, char *s)
3087 {
3088   if (!s) return 0;
3089   else if (p<=0)return 0;
3090   else if (s[p]!='-' && s[p+1]=='-')return 1;
3091   else return 0;
3092 }
3093
3094 char* get_alphabet   ( char *seq, char *alphabet)
3095 /*This function builds an alphabet from the string seq*/
3096 /*Alphabet does not need to be emppty. The total number of unique symbols is returned*/
3097 /*Gaps, as defined in the GAP_LIST are ignored*/
3098     {
3099         int a;
3100         int len;
3101
3102         if ( !alphabet) alphabet=vcalloc ( 200, sizeof (char));
3103
3104         len=strlen (seq);
3105
3106         for ( a=0; a<len; a++)
3107           if ( !is_gap ( seq[a]) && !is_in_set ( seq[a], alphabet+1)){alphabet[(int)++alphabet[0]]=seq[a];}
3108         return alphabet;
3109     }
3110
3111 int array_is_in_set ( char *x, char *set )
3112          {
3113          int len;
3114          int a;
3115
3116          len=strlen (x);
3117          for ( a=0; a< len; a++)
3118              if ( !is_in_set ( x[a], set))return 0;
3119          return 1;
3120          }
3121
3122
3123 int is_in_set ( char r, char *list)
3124         {
3125         static char s[2];
3126
3127         s[0]=r;
3128
3129
3130         if ( strstr ( list,s)!=NULL)return 1;
3131
3132         return 0;
3133         }
3134 char * generate_void ( int x)
3135     {
3136     return generate_string (x, ' ');
3137     }
3138
3139 char * generate_null ( int x)
3140     {
3141     return generate_string ( x, '-');
3142     }
3143
3144 char * generate_string ( int x, char y)
3145     {
3146          int a;
3147     char *string;
3148
3149     string = vcalloc ( x+1, sizeof ( char));
3150
3151     for ( a=0; a< x; a++)
3152         string[a]=y;
3153
3154     string[a]='\0';
3155     return string;
3156     }
3157 char * translate_string (char *string, char *in, char*out)
3158 {
3159   int l1, l2;
3160   int a, b;
3161
3162   l1=strlen(string);
3163   l2=strlen (in);
3164
3165   if ( l2!=strlen (out))
3166     {
3167       fprintf ( stderr, "\nERROR: translate_string function: IN alphabet (%s) must have the same size as OUT alphabet (%s) [FATAL:%s]\n", in, out, PROGRAM);
3168       myexit (EXIT_FAILURE);
3169     }
3170
3171   for ( a=0; a< l1; a++)
3172     {
3173       for ( b=0; b< l2; b++)
3174         {
3175           if ( string[a]==in[b])
3176             {
3177               string[a]=out[b];
3178               b=l2;
3179             }
3180         }
3181     }
3182   return string;
3183 }
3184
3185
3186 int get_longest_string (char **array,int n, int *len, int *index)
3187      {
3188      int a, l;
3189      int max_len=0, max_index=0;
3190
3191      if ( n==0 || array==NULL )return 0;
3192      if ( n==-1)n=read_size_char(array,sizeof (char*));
3193
3194
3195      if (read_size_char ( array,sizeof (char*))<n)
3196          {
3197          fprintf ( stderr, "\nBAD REQUEST: Number of strings=%d, expected number=%d", read_size_char ( array,sizeof (char*)),n);
3198          crash ("[FATAL/util.c]");
3199          }
3200      else
3201          {
3202          max_len=strlen(array[0]);
3203          max_index=0;
3204          for ( a=1; a< n; a++)
3205              {
3206              l=strlen ( array[a]);
3207              if ( l>max_len)
3208                 {
3209                 max_len=l;
3210                 max_index=a;
3211                 }
3212              }
3213          if (index!=NULL)index[0]=max_index;
3214          if (len!=NULL)len[0]=max_len;
3215          }
3216
3217      return max_len;
3218      }
3219
3220 int get_shortest_string (char **array,int n, int *len, int *index)
3221      {
3222      int a, l;
3223      int min_len;
3224
3225      if ( n==0|| array==NULL || read_size_char ( array,sizeof (char*))<n)return 0;
3226      else
3227          {
3228          min_len=strlen(array[0]);
3229
3230          for ( a=1; a< n; a++)
3231              {
3232              l=strlen ( array[a]);
3233              if ( l<min_len)
3234                 {
3235                 min_len=l;
3236
3237                 }
3238              }
3239          if (index!=NULL)index[0]=a;
3240          if (len!=NULL)len[0]=min_len;
3241          }
3242      return min_len;
3243      }
3244
3245 char **pad_string_array ( char **array, int n, int len, char pad)
3246      {
3247      int a, b, l;
3248      for ( a=0; a< n; a++)
3249          {
3250          l= strlen ( array[a]);
3251          for (b=l; b<len; b++)array[a][b]=pad;
3252          array[a][len]='\0';
3253          }
3254      return array;
3255      }
3256 char **right_pad_string_array ( char **array, int n, int len, char pad)
3257      {
3258      return pad_string_array(array, n, len, pad);
3259      }
3260 char **left_pad_string_array ( char **array, int n, int len, char pad)
3261      {
3262      int a, b, c, l;
3263      char *buf;
3264
3265      buf=vcalloc ( len+1, sizeof (char));
3266      for ( a=0; a< n; a++)
3267          {
3268          l= strlen ( array[a]);
3269          for (b=0; b<(len-l); b++)buf[b]=pad;
3270          for (c=0,b=(len-l); b<=len; c++,b++)buf[b]=array[a][c];
3271          sprintf (array[a], "%s", buf);
3272
3273          }
3274      vfree(buf);
3275      return array;
3276      }
3277 char * crop_string (char *string, int start, int end)
3278      {
3279      static char *buf;
3280      /*Extract [start-end[*/
3281
3282      if ( strlen (string)<end)
3283              {
3284              fprintf ( stderr, "\nString=%s, start=%d end=%d", string, start, end);
3285              crash ( "Wrong End in crop String [FATAL]");
3286              }
3287      else
3288          {
3289          buf=vcalloc (strlen (string)+1, sizeof (char));
3290          string[end]='\0';
3291          sprintf ( buf, "%s", string+start);
3292          sprintf ( string,"%s", buf);
3293          vfree (buf);
3294          }
3295
3296      return string;
3297      }
3298
3299 int get_distance2char ( char *x, char*list)
3300     {
3301
3302     char *y;
3303     y=x;
3304     if (x==NULL) return 0;
3305     while (!is_in_set (y[0], list) && y[0]!='\0')y++;
3306     return y-x;
3307     }
3308 /*********************************************************************/
3309 /*                                                                   */
3310 /*                         TIME        FUNCTIONS                     */
3311 /*                                                                   */
3312 /*                                                                   */
3313 /*********************************************************************/
3314 static long ref;
3315 static int child;
3316 static long ticks;
3317 static long milli_sec_conv=1000;
3318
3319
3320 FILE *print_program_information (FILE *fp, char *comment)
3321 {
3322    fprintf ( fp, "# Results Produced with %s (%s)\n",PROGRAM, VERSION);
3323    fprintf ( fp, "# %s is available from %s\n",PROGRAM,URL);
3324    if (comment)fprintf ( fp, "# %s\n", comment);
3325    return fp;
3326 }
3327 FILE* print_cpu_usage (FILE *fp, char *comment)
3328 {
3329   fprintf ( fp, "# %s CPU Usage: %d millisec\n", comment, get_time());
3330   return fp;
3331 }
3332 void print_exit_success_message ()
3333 {
3334
3335   fprintf ( stderr, "\n# TERMINATION STATUS: SUCCESS [PROGRAM: %s]\n", PROGRAM);
3336 }
3337 void print_exit_failure_message ()
3338 {
3339   fprintf ( stderr, "\n# TERMINATION STATUS: FAILURE [PROGRAM: %s]\n", PROGRAM);
3340 }
3341
3342
3343
3344
3345 int get_time ()
3346         {
3347         static long time;
3348         struct tms time_buf[1];
3349         long tms_stime, tms_utime;
3350
3351         if ( child==1)return get_ctime();
3352
3353         if ( ticks==0)ticks = sysconf(_SC_CLK_TCK);
3354         times ( time_buf);
3355
3356         tms_stime=(long)time_buf->tms_stime*milli_sec_conv;
3357         tms_utime=(long)time_buf->tms_utime*milli_sec_conv;
3358
3359
3360
3361
3362         if ( ref==0)
3363                 {
3364                 ref=(tms_stime+tms_utime);
3365                 return 0;
3366                 }
3367         else
3368                 {
3369                 time=(tms_utime+tms_stime)-ref;
3370                 return (int) ((time)/ticks);
3371                 }
3372         }
3373 int get_ctime ()
3374         {
3375         static long time;
3376         struct tms time_buf[1];
3377         long   tms_cutime, tms_cstime;
3378
3379         if ( ticks==0)ticks = sysconf(_SC_CLK_TCK);
3380         times ( time_buf);
3381
3382
3383
3384         tms_cstime=(long)time_buf->tms_cstime*milli_sec_conv;
3385         tms_cutime=(long)time_buf->tms_cutime*milli_sec_conv;
3386
3387         if ( ref==0)
3388                 {
3389                 child=1;
3390                 ref=tms_cstime+tms_cutime;
3391                 return 0;
3392                 }
3393         else
3394                 {
3395                 time=(tms_cutime+tms_cstime)-ref;
3396                 return (int)((time)/ticks);
3397                 }
3398         }
3399 int reset_time()
3400         {
3401         ref=0;
3402         return (int)get_time();
3403         }
3404 int increase_ref_time(int increase)
3405         {
3406         if ( ref==0)get_time();
3407
3408         ref-=(long)ticks*(long)increase;
3409         if (ref==0)ref++;
3410         return (int)ref;
3411         }
3412
3413 /*********************************************************************/
3414 /*                                                                   */
3415 /*                         SYSTEM CALLS                              */
3416 /*                                                                   */
3417 /*                                                                   */
3418 /*********************************************************************/
3419 int evaluate_sys_call_io ( char *out_file, char *com, char *fonc)
3420      {
3421      if ( check_file_exists (out_file))return 1;
3422      else
3423         {
3424         fprintf ( stderr, "\nCommand\n%s\nFailed to produce File\n%s\n", com, out_file);
3425         return 0;
3426         }
3427      }
3428 void HERE (char *string, ...)
3429 {
3430   va_list ap;
3431
3432   va_start (ap, string);
3433   fprintf ( stderr, "HERE: ");
3434   vfprintf (stderr, string, ap);
3435   fprintf ( stderr, "\n");
3436   va_end (ap);
3437
3438 }
3439 void printf_exit  (int exit_code, FILE *fp, char *string, ...)
3440 {
3441
3442
3443   va_list ap;
3444
3445   va_start (ap, string);
3446   vfprintf (fp, string, ap);
3447   va_end (ap);
3448   myexit (exit_code);
3449 }
3450
3451
3452 int fprintf_fork  (FILE *fp, char *string, ...)
3453 {
3454   va_list ap;
3455   static char *openF;
3456   static char *closeF;
3457
3458   char *pid_file;
3459   FILE *flag;
3460   struct flock fl;
3461   int fd,a;
3462
3463   char buf[100000];
3464
3465   if (!openF)
3466     {
3467       openF=vcalloc (100, sizeof (char));
3468       sprintf (openF, "cedric1");
3469       closeF=vcalloc (100, sizeof (char));
3470       sprintf (closeF, "cedric2");
3471
3472       //openF =vtmpnam (NULL);
3473       //closeF=vtmpnam (NULL);
3474       vfclose(vfopen (openF,"w"));
3475     }
3476   while ((rename (openF,closeF))==-1);
3477
3478   va_start (ap, string);
3479   vsprintf (buf, string, ap);
3480   va_end (ap);
3481   fprintf ( fp, "%s", buf);
3482   fflush (fp);
3483   rename (closeF, openF);
3484
3485   return 0;
3486 }
3487 int fprintf_fork2  (FILE *fp, char *string, ...)
3488 {
3489   va_list ap;
3490   static char *openF;
3491   static struct flock fl;
3492   char buf[100000];
3493   int fd;
3494   va_start (ap, string);
3495   vsprintf (buf, string, ap);
3496   va_end (ap);
3497
3498   fprintf ( fp, "%s", buf);
3499   fflush (fp);
3500   return 0;
3501 }
3502
3503 int printf_file  (char *file,char *mode, char *string,...)
3504 {
3505   FILE *fp;
3506   va_list ap;
3507
3508   if (!(fp=vfopen (file, mode)))return 0;
3509   va_start (ap, string);
3510   vfprintf (fp, string, ap);
3511   va_end (ap);
3512   vfclose (fp);
3513   return 1;
3514   }
3515 int printf_system_direct  (char *string, ...)
3516 {
3517   char buf[10000];
3518
3519   va_list ap;
3520
3521   va_start (ap, string);
3522   vsprintf (buf, string, ap);
3523   va_end (ap);
3524   return safe_system (buf);
3525 }
3526
3527 int printf_system  (char *string, ...)
3528 {
3529   char buf[10000];
3530
3531   va_list ap;
3532
3533   va_start (ap, string);
3534   vsprintf (buf, string, ap);
3535   va_end (ap);
3536   return my_system (buf);
3537 }
3538
3539 int my_system_cl (int argc, char *argv[])
3540 {
3541   int a,l;
3542   char *command;
3543
3544   for ( a=0, l=0; a< argc; a++)l+=(strlen(argv[a])+2);
3545   command=vcalloc (l+1, sizeof(char));
3546   for ( a=0; a< argc; a++)
3547     {
3548       command=strcat (command, argv[a]);
3549       command=strcat (command, " ");
3550     }
3551   a=my_system ( command);
3552   vfree (command);
3553   return a;
3554 }
3555
3556 int my_system ( char *command0)
3557 {
3558   static char ***unpacked_list;
3559   static int n_unpacked;
3560   static int proxy_set;
3561   static int email_set;
3562   int email=0, proxy=0, update_env=0;
3563
3564   //Set the net and E-mail status
3565   //if ( strstr (command0, "wget"))proxy=1;
3566   //if ( strstr (command0, "curl"))proxy=1;
3567
3568
3569
3570   if ( strstr (command0, "extract_from_pdb"))proxy=1;
3571   if ( strstr (command0, "tc_generic_method"))proxy=1;
3572   if ( strstr (command0, "install"))proxy=1;
3573   if ( strstr (command0, "tc_generic_method"))email=1;
3574
3575
3576
3577   if (!unpacked_list)
3578     {
3579       unpacked_list=declare_arrayN(3, sizeof (char), 3, 200,300);
3580     }
3581
3582   if ( getenv ("DEBUG_PERL"))return safe_system (command0);
3583   else
3584     {
3585       char **list;
3586       int is_command;
3587       int a, c=0;
3588       char *command1;
3589       char *command2;
3590       int return_val;
3591
3592       command1=vcalloc ( 3*strlen (command0)+1, sizeof (char));
3593       command2=vcalloc ( 100000, sizeof (char));
3594       sprintf ( command1, "%s", command0);
3595
3596       command1=substitute (command1, "|", " | ");
3597       command1=substitute (command1, ";", " ; ");
3598
3599       list=string2list (command1);
3600       if ( !list) return EXIT_SUCCESS;
3601       is_command=1;
3602
3603       //Identify T-Coffee self threads and install threads
3604       if ( strstr (list[1], "t_coffee"))update_env=1;
3605       else if ( strstr (list[1], "install.pl"))check_internet_connection (IS_FATAL);
3606
3607
3608       for ( a=1; a< atoi(list[0]); a++)
3609         {
3610           if ( is_command)
3611             {
3612               if ( strstr ( list[a], "unpack_"))
3613                 {
3614                   unpack_all_perl_script (list[a]+strlen ("unpack_"));
3615                   myexit (EXIT_SUCCESS);
3616                 }
3617               else if ((c=name_is_in_list (list[a], unpacked_list[0], n_unpacked, 100))!=-1);
3618               else
3619                 {
3620                   n_unpacked=unpack_perl_script (list[a], unpacked_list, n_unpacked);c=n_unpacked-1;
3621                 }
3622               //if non unpacked script check pg is installed:
3623
3624               if ( strm (unpacked_list[2][c], "shell"))
3625                 {
3626                   check_program_is_installed (list[a], NULL, NULL, NULL, INSTALL_OR_DIE);
3627                 }
3628               strcat (command2, ((c!=-1)?unpacked_list[1][c]:list[a]));
3629               strcat (command2, " ");
3630               is_command=0;
3631
3632             }
3633           else
3634             {
3635               strcat (command2, list[a]);
3636               strcat (command2, " ");
3637               if ( strm (list[a], ",") ||strm (list[a], "|")) is_command=1;
3638             }
3639         }
3640
3641       free_char (list,-1);
3642       vfree ( command1);
3643       command2=substitute ( command2, "//", "/");
3644
3645       return_val=safe_system (command2);
3646
3647       //INTERCEPT POTENTIAL NETWORK CALLS
3648
3649       if (return_val!=EXIT_SUCCESS && proxy==1 && proxy_set==0)
3650         {
3651           if (simple_check_internet_connection (NULL))proxy_set=1;
3652           else
3653             {
3654               check_internet_connection (IS_NOT_FATAL);proxy_set=1;
3655               return_val=safe_system (command2);
3656             }
3657         }
3658       //Intercept potential Missing E-mails
3659       if (return_val!=EXIT_SUCCESS && email==1 && email_set==0)
3660         {
3661           Email(INPUT, RESET);email_set=1;
3662           return_val=safe_system (command2);
3663         }
3664
3665       //update the environement that may have been modified by a thread
3666
3667       if (update_env)get_t_coffee_environement (NULL);
3668
3669       vfree ( command2);
3670       return return_val;
3671     }
3672 }
3673 int safe_system (const char * com)
3674 {
3675     pid_t pid;
3676     int   status;
3677     if (com == NULL)
3678         return (1);
3679
3680     if ((pid = fork ()) < 0)
3681         return (-1);
3682
3683     if (pid == 0) {
3684
3685         char * argv [4];
3686
3687         argv [0] = "sh";
3688         argv [1] = "-c";
3689         argv [2] =(char*) com;
3690         argv [3] = 0;
3691         execvp ("/bin/sh", argv);
3692     }
3693     else
3694       {
3695         set_pid(pid);
3696       }
3697
3698
3699     while (1) {
3700
3701         if (vwaitpid (pid, &status, 0) == -1)
3702           {
3703             if (errno != EINTR)
3704               return (EXIT_FAILURE);
3705           }
3706         else
3707           {
3708             return (status);
3709           }
3710     }
3711 }
3712 static int **pidtable;
3713 static int pidtable_s;
3714 pid_t **declare_pidtable ()
3715 {
3716   int a;
3717   pidtable_s=MAX_N_PID;
3718
3719   pidtable=vcalloc (pidtable_s, sizeof (pid_t*));
3720   for (a=0; a< pidtable_s; a++)
3721     {
3722       pidtable[a]=vcalloc (2, sizeof (pid_t));
3723     }
3724   return pidtable;
3725 }
3726 pid_t set_pid (pid_t p)
3727 {
3728   int cpid;
3729
3730   if (!pidtable)declare_pidtable();
3731   if ( p<=0) return;
3732   else if ( p>=MAX_N_PID)printf_exit (EXIT_FAILURE,stderr,"ERROR PID=%d superior to MAX_N_PID=%d [FATAL]",p, MAX_N_PID);
3733   pidtable[(int)p][0]=getpid();
3734   pidtable[(int)p][1]=1;
3735 }
3736 pid_t vfork ()
3737 {
3738   pid_t p;
3739   static int attempt;
3740
3741   if ( attempt==1000) printf_exit (EXIT_FAILURE, stderr,"\nERROR: Could not fork processes. Run again with -multi_core=no\n");
3742
3743
3744   p=fork();
3745   if (p==-1)
3746     {
3747       attempt++;
3748       wait(-1);
3749       return vfork();
3750     }
3751   else
3752     {
3753       attempt=0;
3754       return p;
3755     }
3756 }
3757 int    vwait_npid (int sub, int max, int min)
3758 {
3759   if (max==0)
3760     {
3761       while (sub>0)
3762         {
3763           vwait (NULL);
3764           sub--;
3765         }
3766     }
3767   else if ( sub>=max)
3768     {
3769       while (sub>=min)
3770         {
3771           vwait (NULL);
3772           sub--;
3773         }
3774     }
3775   else{;}
3776   return sub;
3777 }
3778
3779 pid_t  vwaitpid (pid_t p, int *status, int options)
3780 {
3781   pid_t p2;
3782
3783   p=waitpid (p, status, options);
3784
3785   if (pidtable)pidtable[(int)p][0]=pidtable[(int)p][1]=0;
3786   return p;
3787 }
3788 pid_t vwait (pid_t *p)
3789 {
3790   pid_t p2;
3791
3792   p2=wait (p);
3793
3794   if (pidtable)pidtable[(int)p2][0]=pidtable[(int)p2][1]=0;
3795   return p2;
3796 }
3797 int kill_child_pid()
3798 {
3799   int n;
3800
3801   if ( !pidtable)return 0;
3802   else
3803     {
3804       int a;
3805       pid_t cpid;
3806       cpid=getpid();
3807       for (a=0; a<pidtable_s; a++)
3808         {
3809           if (pidtable[a][1] && pidtable[a][0]==cpid )
3810             {
3811               pidtable[a][1]=pidtable[a][0]=0;
3812               kill((pid_t)a, SIGTERM);n++;
3813             }
3814         }
3815     }
3816   return n;
3817 }
3818
3819
3820
3821
3822 void unpack_all_perl_script (char *script)
3823 {
3824   int a=0;
3825   FILE *fp;
3826   char command[1000];
3827
3828   while ( !strm(PerlScriptName[a], "EndList"))
3829     {
3830       if (script==NULL || strm (script, "all") ||strm (script, "perl_script")|| strm (script,PerlScriptName[a]))
3831         {
3832           fprintf ( stderr, "\nUnpack Script %s\n", PerlScriptName[a]);
3833           fp=vfopen ( PerlScriptName[a], "w");
3834           fprintf (fp, "%s\n%s\n", PERL_HEADER,PerlScriptFile[a]);
3835           vfclose (fp);
3836           sprintf ( command, "chmod u+x %s",PerlScriptName[a] );
3837           safe_system (command);
3838         }
3839       a++;
3840     }
3841   myexit (EXIT_SUCCESS);
3842 }
3843 int satoi (char *);
3844 int satoi (char *c)
3845 {
3846   if ( !c)return 0;
3847   else return atoi (c);
3848 }
3849
3850 int unpack_perl_script (char *name, char ***unpacked, int n)
3851 {
3852   int a=0;
3853
3854   set_file2remove_extension(".pl", SET);
3855
3856   if ( name==NULL) unpack_all_perl_script (NULL);
3857   while ( !strm(PerlScriptName[a], "EndList") && !strm ( name, PerlScriptName[a]))
3858     {
3859       a++;
3860     }
3861
3862   if ( strm(PerlScriptName[a], "EndList")|| (check_file_exists (name) && isexec(name) && satoi(getenv("RUN_LOCAL_SCRIPT"))) || (pg_is_installed (name) &&(satoi(getenv("RUN_LOCAL_SCRIPT"))) ))
3863        {
3864          if ( check_file_exists (name) && ! strstr (name, "/") && isexec(name) && satoi(getenv ("RUN_LOCAL_SCRIPT")))
3865            {
3866              add_warning ( stderr, "\nWARNING: Local Copy of %s detected [Debug mode]. DELETE this file\n", name);
3867              sprintf ( unpacked[0][n], "%s", name);
3868              sprintf ( unpacked[1][n], "./%s", name);
3869              sprintf ( unpacked[2][n], "shell");
3870            }
3871          else if ( pg_is_installed (name) &&satoi(getenv("RUN_LOCAL_SCRIPT")) )
3872            {
3873              add_warning ( stderr, "\nWARNING: Running external copy of %s [Debug mode].\n", name);
3874              sprintf ( unpacked[0][n], "%s", name);
3875              sprintf ( unpacked[1][n], "%s", name);
3876              sprintf ( unpacked[2][n], "local");
3877            }
3878          else if ( strstr (name, "t_coffee"))
3879            {
3880              /*Cygwin pb: t_coffee cannot call itself: making a local copy for recursion*/
3881              char buf[1000];
3882              char *b2;
3883              b2=vtmpnam (NULL);
3884              sprintf (buf, "cp `which %s` %s; chmod u+x %s",name, b2, b2);
3885              safe_system (buf);
3886              sprintf ( unpacked[0][n], "%s", name);
3887              sprintf ( unpacked[1][n], "%s", b2);
3888              sprintf ( unpacked[2][n], "shell");
3889
3890            }
3891          else
3892            {
3893              sprintf ( unpacked[0][n], "%s", name);
3894              sprintf ( unpacked[1][n], "%s", name);
3895              sprintf ( unpacked[2][n], "shell");
3896            }
3897        }
3898   else
3899     {
3900       FILE *fp;
3901       char command[1000];
3902       sprintf ( unpacked[0][n], "%s", name);
3903       sprintf ( unpacked[1][n], "%s", vtmpnam(NULL));
3904       sprintf ( unpacked[2][n], "unpack");
3905       fp=vfopen (unpacked[1][n], "w");
3906       fprintf (fp, "%s\n%s\n", PERL_HEADER,PerlScriptFile[a]);
3907       vfclose (fp);
3908       sprintf ( command, "chmod u+x %s", unpacked[1][n]);
3909       safe_system (command);
3910     }
3911   set_file2remove_extension(".pl", UNSET);
3912   return ++n;
3913 }
3914
3915 /*********************************************************************/
3916 /*                                                                   */
3917 /*                         IO FUNCTIONS                              */
3918 /*                                                                   */
3919 /*                                                                   */
3920 /*********************************************************************/
3921 char *program_name;
3922 char *error_file;
3923 static int add_file2error_report ( char *type, char *format, char *name);
3924
3925 char *in_cl;
3926 FILE * print_command_line (FILE *fp )
3927 {
3928   fprintf ( fp, "# Command Line: %s [PROGRAM:%s]\n",in_cl, PROGRAM);
3929   return fp;
3930 }
3931 int check_dir_getenv ( char *string);
3932 char *get_plugins_4_tcoffee ( char *mode)
3933 {
3934   static char *plugins_4_tcoffee;
3935
3936   if ( !plugins_4_tcoffee)plugins_4_tcoffee=vcalloc ( 1000, sizeof (char));
3937
3938   if (mode)plugins_4_tcoffee[0]='\0';
3939
3940
3941   if ( plugins_4_tcoffee[0])return plugins_4_tcoffee;
3942   else if ( mode) sprintf (plugins_4_tcoffee, "%s", mode);
3943   else if ( check_dir_getenv ("PLUGINS_4_TCOFFEE"))
3944     {
3945       sprintf (plugins_4_tcoffee, "%s", getenv ("PLUGINS_4_TCOFFEE"));
3946     }
3947   else
3948     {
3949       sprintf (plugins_4_tcoffee, "%s/plugins/%s/", get_dir_4_tcoffee(), get_os());
3950     }
3951   //if (!isdir (plugins_4_tcoffee)){printf_exit (EXIT_FAILURE, stderr, "ERROR: The declared plugin directory for T-Coffee is invalid: [%s] [FATAL:%s %s]", plugins_4_tcoffee, PROGRAM, VERSION); }
3952   return plugins_4_tcoffee;
3953 }
3954 char *get_home_4_tcoffee ()
3955 {
3956   static char *home_4_tcoffee;
3957   static char home[1000];
3958
3959
3960   if ( !home_4_tcoffee)
3961     home_4_tcoffee=vcalloc ( 1000, sizeof (char));
3962
3963   if ( home_4_tcoffee[0])return home_4_tcoffee;
3964   else if ( check_dir_getenv ("HOME_4_TCOFFEE"))
3965     {
3966       sprintf (home_4_tcoffee, "%s", getenv ("HOME_4_TCOFFEE"));
3967     }
3968   else if ( check_dir_getenv ("HOME"))
3969     {
3970       sprintf (home_4_tcoffee, "%s", getenv ("HOME"));
3971       sprintf (home, "%s", home_4_tcoffee);
3972     }
3973   else if ( check_dir_getenv ("TMP"))
3974     {
3975       sprintf (home_4_tcoffee, "%s", getenv ("TMP"));
3976     }
3977   else if (check_dir_getenv ("TEMP"))
3978     {
3979       sprintf (home_4_tcoffee, "%s", getenv ("TEMP"));
3980     }
3981   else
3982     {
3983       printf_exit (EXIT_FAILURE, stderr, "ERROR: Could not set a HOME directory.\nSet any of the following environement variables to some suitable location: HOME, HOME_4_TCOFFEE, TMP or TEMP [FATAL:%s]\n", PROGRAM);
3984
3985     }
3986
3987
3988   if ( !strm( home, home_4_tcoffee))
3989     add_warning (stderr, "WARNING: The Home directory was set to %s", home_4_tcoffee);
3990
3991   return home_4_tcoffee;
3992 }
3993 int get_nproc ()
3994 {
3995   static int nproc;
3996
3997   if (nproc);
3998   else if ( get_int_variable ("n_core"))
3999     {
4000       nproc=get_int_variable ("n_core");
4001     }
4002   else if ( getenv ("NUMBER_OF_PROCESSORS_4_TCOFFEE"))
4003     {
4004       nproc=atoi(getenv ("NUMBER_OF_PROCESSORS_4_TCOFFEE"));
4005     }
4006   else if ( getenv ("NUMBER_OF_PROCESSORS"))
4007     {
4008       nproc=atoi(getenv ("NUMBER_OF_PROCESSORS"));
4009     }
4010   else if ( check_file_exists ("/proc/cpuinfo"))
4011     {
4012       char *tmp;
4013       char string[100];
4014       tmp=vtmpnam (NULL);
4015       printf_system ("grep \"^processor\" /proc/cpuinfo | wc -l>%s", tmp);
4016       sprintf ( string, "%s", file2string (tmp));
4017       chomp (string);
4018       nproc=atoi (string);
4019     }
4020   else
4021     nproc=1;
4022
4023   return nproc;
4024 }
4025 char * get_os()
4026 {
4027   static char os[100];
4028   char *file;
4029
4030   if ( os[0])return os;
4031   else
4032     {
4033       char *command;
4034       char *s;
4035
4036       command=vcalloc (100, sizeof (char));
4037       file=tmpnam (NULL);
4038       sprintf ( command, "uname > %s", file);
4039       safe_system (command);
4040       s=file2string (file);
4041       lower_string (s);
4042
4043       if (strstr (s, "cygwin"))sprintf ( os, "windows");
4044       else if ( strstr (s, "linux"))sprintf ( os, "linux");
4045       else if ( strstr (s, "osx"))sprintf ( os, "macosx");
4046       else if ( strstr (s, "darwin"))sprintf ( os, "macosx");
4047       else sprintf (os, "%s", s);
4048       vfree (s);
4049       vfree (command);
4050       vremove (file);
4051     }
4052   return os;
4053 }
4054
4055 char *file_putenv (char *file)
4056 {
4057   if (!file) return NULL;
4058   else if ( !check_file_exists (file))return NULL;
4059   else
4060     {
4061       char ***list;
4062       int n=0;
4063       list=file2list (file, "\n");
4064       while (list[n])
4065         {
4066           if ( list[n][1][0]!='#')cputenv ("%s",list[n][1]);
4067           n++;
4068         }
4069       free_arrayN ((void ***)list, 3);
4070     }
4071
4072   return NULL;
4073 }
4074 int cputenv   (char * string, ...)
4075 {
4076   va_list ap;
4077   char *file;
4078   char *s;
4079   FILE *fp;
4080   if (!string)return 0;
4081   file=vtmpnam (NULL);
4082   va_start (ap, string);
4083   fp=vfopen (file, "w");
4084   vfprintf (fp, string, ap);
4085   vfclose (fp);
4086   va_end (ap);
4087   s=file2string (file);
4088   if (!s) return 0;
4089   putenv (s);
4090   return 1;
4091 }
4092
4093 int check_dir_getenv ( char *string)
4094 {
4095   char *p;
4096
4097
4098
4099   p=getenv ( string);
4100   if ( !p) return 0;
4101   if ( !p || access (p, F_OK)==-1 || access (p, W_OK)==-1 || access(p, R_OK)==-1 || access (p, X_OK)==-1)return 0;
4102
4103   return 1;
4104 }
4105
4106 char *get_dir_4_tcoffee()
4107 {
4108   static char dir_4_tcoffee[1000];
4109   if (dir_4_tcoffee[0])return dir_4_tcoffee;
4110   else
4111     {
4112       if ( getenv ("DIR_4_TCOFFEE"))sprintf (dir_4_tcoffee, "%s", getenv("DIR_4_TCOFFEE"));
4113       else sprintf ( dir_4_tcoffee, "%s/.t_coffee",get_home_4_tcoffee());
4114       sprintf ( DIR_4_TCOFFEE, "%s", dir_4_tcoffee);
4115       my_mkdir (dir_4_tcoffee);
4116     }
4117   return dir_4_tcoffee;
4118 }
4119 char *get_tmp_4_tcoffee ()
4120 {
4121   static char tmp_4_tcoffee [1000];
4122
4123   if ( tmp_4_tcoffee[0])return tmp_4_tcoffee;
4124   else
4125     {
4126
4127       if ( getenv ("TMP_4_TCOFFEE"))sprintf (tmp_4_tcoffee, "%s", getenv("TMP_4_TCOFFEE"));
4128       else
4129         {
4130           char command [1000];
4131
4132           if ( strm (get_os(), "windows"))
4133             {
4134               sprintf ( tmp_4_tcoffee, ".TCtmp");
4135             }
4136           else
4137             {
4138             sprintf ( tmp_4_tcoffee, "%s/tmp", get_dir_4_tcoffee());
4139             }
4140         }
4141
4142       sprintf ( TMP_4_TCOFFEE, "%s", tmp_4_tcoffee);
4143       my_mkdir (tmp_4_tcoffee);
4144     }
4145
4146   return tmp_4_tcoffee;
4147 }
4148 char *get_cache_4_tcoffee ()
4149 {
4150
4151   static char cache_4_tcoffee [1000];
4152   if ( cache_4_tcoffee[0])return cache_4_tcoffee;
4153   else
4154     {
4155       if ( getenv ("CACHE_4_TCOFFEE"))sprintf (cache_4_tcoffee, "%s", getenv("CACHE_4_TCOFFEE"));
4156       else sprintf ( cache_4_tcoffee, "%s/cache/", get_dir_4_tcoffee());
4157       sprintf ( CACHE_4_TCOFFEE, "%s", cache_4_tcoffee);
4158       my_mkdir(cache_4_tcoffee); /*Do not use mkdir: not yet initialized*/
4159     }
4160   return cache_4_tcoffee;
4161 }
4162 char *get_mcoffee_4_tcoffee ()
4163 {
4164   static char mcoffee_4_tcoffee [1000];
4165   if ( mcoffee_4_tcoffee[0])return mcoffee_4_tcoffee;
4166   else
4167     {
4168       if ( getenv ("MCOFFEE_4_TCOFFEE"))sprintf (mcoffee_4_tcoffee, "%s", getenv("MCOFFEE_4_TCOFFEE"));
4169       else sprintf ( mcoffee_4_tcoffee, "%s/mcoffee/", get_dir_4_tcoffee());
4170       my_mkdir (mcoffee_4_tcoffee);
4171     }
4172   return mcoffee_4_tcoffee;
4173 }
4174 char *get_methods_4_tcoffee ()
4175 {
4176   static char methods_4_tcoffee [1000];
4177   if ( methods_4_tcoffee[0])return methods_4_tcoffee;
4178   else
4179     {
4180       if ( getenv ("METHODS_4_TCOFFEE"))sprintf (methods_4_tcoffee, "%s", getenv("METHODS_4_TCOFFEE"));
4181       else sprintf ( methods_4_tcoffee, "%s/methods/", get_dir_4_tcoffee());
4182       my_mkdir(methods_4_tcoffee);
4183     }
4184   return methods_4_tcoffee;
4185 }
4186
4187 char ** standard_initialisation ( char **in_argv, int *in_argc)
4188 {
4189   return standard_initialisation_end   (standard_initialisation_start (in_argv, in_argc), in_argc);
4190 }
4191 int getpid_ref()
4192 {
4193   static int pid;
4194   if (!pid)pid=getpid();
4195   return pid;
4196 }
4197
4198 char ** standard_initialisation_start ( char **in_argv, int *in_argc)
4199       {
4200         static int standard_initialisation_done;
4201         char **out_argv;
4202         FILE *fp, *fp2;
4203         int a,c, stdi;
4204         //Break the command line, intercept the pipe, prepare the exit
4205
4206         getpid_ref();
4207         if ( in_argv==NULL)
4208           {
4209             standard_initialisation_done=0;
4210             return NULL;
4211           }
4212
4213         else if ( standard_initialisation_done==1)
4214           {
4215             return in_argv;
4216           }
4217         else standard_initialisation_done=1;
4218         /*1 Check for the cache and tmp directories*/
4219
4220         get_dir_4_tcoffee();
4221         get_tmp_4_tcoffee();
4222
4223         for (c=0,a=0; a<in_argc[0]; a++)c+=strlen (in_argv[a])+3;
4224         in_cl=vcalloc ( c, sizeof (char));
4225         for (a=0; a<in_argc[0]; a++)
4226           {
4227             strcat ( in_cl, in_argv[a]);
4228             strcat (in_cl, " ");
4229           }
4230
4231 /*Standard exit*/
4232         global_exit_signal=EXIT_SUCCESS;
4233         atexit (main_exit);
4234         atexit (output_warning_list);
4235         signal (SIGTERM, main_exit);
4236         signal (SIGINT, main_exit);
4237
4238         signal (SIGABRT ,error_exit);
4239         signal (SIGFPE  ,error_exit);
4240         signal (SIGILL  ,error_exit);
4241         signal (SIGSEGV ,error_exit);
4242
4243
4244 /*Process the parameters*/
4245
4246       program_name=vcalloc ( strlen (in_argv[0])+strlen (PROGRAM)+1, sizeof (char));
4247       if (in_argv)
4248           {
4249             sprintf ( program_name, "%s", in_argv[0]);
4250             out_argv=break_list ( in_argv, in_argc, "=;, \n");
4251           }
4252       else
4253         {
4254          sprintf ( program_name, "%s",PROGRAM);
4255         }
4256
4257       if ( name_is_in_list ( "-no_error_report", out_argv, in_argc[0], 100)!=-1)
4258         no_error_report=1;
4259
4260 /*Error Log File: contains every input file*/
4261       error_file=vtmpnam(NULL);
4262       fp=vfopen (error_file, "w");
4263       fprintf ( fp, "\n######### RUN_REPORT      START  ######");
4264       fprintf ( fp, "\n######### PROGRAM_VERSION START  ######");
4265       fprintf ( fp, "\n          %s, %s", PROGRAM, VERSION);
4266       fprintf ( fp, "\n######### PROGRAM_VERSION END    ######");
4267
4268       fprintf ( fp, "\n######### COMMAND_LINE    START  ######");
4269       fprintf ( fp, "\n");
4270       for (a=0; a<in_argc[0]; a++)fprintf (fp, "%s ", out_argv[a]);
4271       fprintf ( fp, "\n######### COMMAND_LINE    END    ######\n");
4272
4273       for ( a=1; a<in_argc[0]; a++)
4274         {
4275           if ( check_file_exists ( out_argv[a]) && !filename_is_special (out_argv[a]))
4276             {
4277
4278               fprintf ( fp, "\n######### FILE: %s        START  ######\n", out_argv[a]);
4279               fprintf ( fp, "\n");
4280               fp2=vfopen ( out_argv[a], "r");
4281               while ((c=fgetc(fp2))!=EOF)fprintf (fp, "%c", c);
4282               vfclose (fp2);
4283               fprintf ( fp, "\n######### FILE: %s        END    ######\n", out_argv[1]);
4284             }
4285         }
4286       fprintf ( fp, "\n######### RUN_REPORT END         ######\n");
4287       vfclose (fp);
4288
4289 /*Initialise the time*/
4290       get_time();
4291
4292
4293       for (a=1, stdi=0; a<in_argc[0]; a++)
4294         {
4295           if ( (strm ( out_argv[a], "stdin") || strm (out_argv[a], "STDIN")) && stdi==0)
4296             {
4297               char *file;
4298               FILE *fp, *fp_stdi;
4299
4300               if ( stdi==0)
4301                 {
4302                   file=vtmpnam (NULL);
4303                   vfree (out_argv[a]);
4304                   out_argv[a]=vcalloc ( strlen (file)+1, sizeof (char));
4305                   sprintf (out_argv[a], "%s", file);
4306                   fp_stdi=vfopen ("stdin", "r");
4307                   fp=vfopen (file, "w");
4308                   while ( (c=fgetc (fp_stdi))!=EOF)fprintf (fp, "%c", c);
4309                   vfclose (fp);
4310                   vfclose (fp_stdi);
4311                   stdi=1;
4312                 }
4313               else if ( stdi ==1)
4314                 {
4315                   printf_exit ( EXIT_FAILURE,stderr, "ERROR: Attempt to get more than one file via pipe [FATAL:%s]\n", PROGRAM);
4316                 }
4317             }
4318         }
4319
4320       return out_argv;
4321       }
4322 char ** standard_initialisation_end ( char **in_argv, int *in_argc)
4323 {
4324   //All the initialization that
4325   static int standard_initialisation_end_done;
4326   //Break the command line, intercept the pipe, prepare the exit
4327   if ( in_argv==NULL)
4328     {
4329       standard_initialisation_end_done=0;
4330       return NULL;
4331     }
4332   else if ( standard_initialisation_end_done==1)
4333     {
4334       return in_argv;
4335     }
4336   else standard_initialisation_end_done=1;
4337
4338   //set the environement variables (proxy and E-mail if possible)
4339   get_t_coffee_environement (NULL);
4340
4341
4342   //clean the cache if needed
4343   //printf_system ("clean_cache.pl -dir=%s/cache -size=%d",get_dir_4_tcoffee (), CACHE_MAX_SIZE);
4344   return in_argv;
4345 }
4346
4347
4348 int add_file2error_report ( char *type, char *format, char *name)
4349   {
4350   FILE *fp1, *fp2;
4351   int c;
4352
4353
4354
4355   fp1=vfopen ( error_file, "a");
4356   fp2=vfopen (name, "r");
4357   fprintf ( fp1, "######### OUTFILE: %s        START  ######\n", name);
4358   while ( (c=fgetc(fp2))!=EOF)fprintf (fp1, "%c", c);
4359   c=ungetc(c,fp2);c=fgetc(fp2); if ( c!='\n')fprintf (fp1, "\n");
4360   fprintf ( fp1, "######### OUTFILE: %s        END    ######\n", name);
4361   vfclose ( fp1);vfclose(fp2);
4362   return 1;
4363
4364 }
4365
4366
4367 void myexit (int signal)
4368 {
4369   global_exit_signal=signal;
4370   if ( signal == EXIT_FAILURE)
4371     {
4372       print_exit_failure_message();
4373     }
4374   else if ( signal==EXIT_SUCCESS)
4375     {
4376       /*print_exit_success_message();*/;
4377     }
4378   exit (global_exit_signal);
4379 }
4380
4381
4382 static int n_warning;
4383 static char **warning_list;
4384
4385 FILE *fatal_exit (FILE *fp,int exit_mode, char *string, ...)
4386 {
4387   va_list ap;
4388   va_start (ap, string);
4389   vfprintf (fp, string, ap);
4390   va_end (ap);
4391   myexit (exit_mode);
4392   return fp;
4393 }
4394 static int warning_mode;
4395 int set_warning_mode ( int mode)
4396 {
4397   warning_mode=mode;
4398   return mode;
4399 }
4400 FILE *add_warning (FILE *fp, char *string, ...)
4401 {
4402   char buf[10000];
4403   static int max_n_warning=100;
4404   va_list ap;
4405
4406
4407   if ( warning_mode==NO || getenv("NO_WARNING_4_TCOFFEE"))return fp;
4408
4409
4410   if ( !warning_list)warning_list=vcalloc (max_n_warning, sizeof (char*));
4411   else if ( n_warning>=max_n_warning)
4412     {
4413       max_n_warning+=100;
4414       warning_list=vrealloc ( warning_list,sizeof (char*)*max_n_warning);
4415     }
4416
4417   va_start (ap, string);
4418   if (fp)
4419     {
4420       fprintf ( fp, "\n");
4421       vfprintf (fp, string, ap);
4422     }
4423   va_end(ap);
4424
4425   va_start (ap, string);
4426   vsprintf (buf, string, ap);
4427   va_end (ap);
4428
4429   warning_list[n_warning]=vcalloc (strlen (buf)+1, sizeof (char));
4430   sprintf ( warning_list[n_warning], "%s", buf);
4431   n_warning++;
4432
4433   return fp;
4434 }
4435 void output_warning_list()
4436 {
4437   int a;
4438   if ( n_warning==0){;}
4439   else
4440     {
4441       fprintf ( stderr, "\nWARNING RECAPITULATION: %d Warning%c [PROGRAM: %s]\n", n_warning, (n_warning>1)?'s':' ', PROGRAM);
4442       for (a=0; a< n_warning; a++)
4443         {
4444           fprintf (stderr, "**WARNING: %3d** %s\n",a+1, warning_list[a]);
4445         }
4446     }
4447 }
4448
4449
4450 int   count_n_res_in_array  (char *array, int len)
4451       {
4452         return count_n_symbol_in_array(array, "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ", len);
4453       }
4454 int   count_n_gap_in_array  (char *array, int len)
4455       {
4456         int l;
4457         if ( len<=0 ||len>strlen(array) )l=strlen(array);
4458         else l=len;
4459
4460         return l- count_n_res_in_array (array,len);
4461       }
4462 int count_n_symbol_in_array ( char *array, char *array_list, int len)
4463       {
4464         int a=0, t=0;
4465         int l;
4466
4467         if ( len<=0 ||len>strlen(array) )l=strlen(array);
4468         else l=len;
4469
4470         for ( a=0; a< l; a++)t+=is_in_set (array[a], array_list);
4471         return t;
4472       }
4473
4474 char* count_strings_in_file ( char *in, char *out)
4475 {
4476   FILE *fp;
4477   int n,c;
4478   char **list, **result;
4479
4480
4481   if (!out) out=vtmpnam (NULL);
4482   list=declare_char (count_n_line_in_file(in)+1, measure_longest_line_in_file (in)+1);
4483
4484   n=0;
4485   fp=vfopen (in, "r");
4486   while ((c=fgetc (fp))!=EOF)
4487     {
4488       ungetc (c, fp);
4489       fscanf (fp, "%s\n",list[n]);
4490       n++;
4491     }
4492   vfclose (fp);
4493
4494   result=count_strings (list, n);
4495
4496
4497   n=0;
4498   fp=vfopen (out, "w");
4499   while (result[n])fprintf ( fp,"%s\n", result[n++]);
4500   vfclose (fp);
4501
4502   free_char (list, -1);
4503   free_char (result, -1);
4504
4505   return out;
4506 }
4507
4508 int ** count_int_strings (int **array, int len, int s)
4509 {
4510   int **result;
4511   int a,n;
4512
4513   sort_list_int (array,s, s,0, len-1);
4514   result=vcalloc (len, sizeof (int*));
4515   for (n=-1,a=0; a<len; a++)
4516     {
4517       if (n==-1 || memcmp (array[a], result[n], sizeof (int)*s)!=0)
4518         {
4519           n++;
4520           result[n]=vcalloc (s+1, sizeof (int));
4521           memcpy (result[n],array[a], sizeof (int)*s);
4522           fprintf ( stdout, "*");
4523         }
4524       result[n][s]++;
4525     }
4526   sort_int (result, s+1, s, 0, n);
4527   result=vrealloc ( result, sizeof (int*)*n);
4528   return result;
4529 }
4530
4531 char** count_strings ( char **array, int len)
4532 {
4533   int a, c, n;
4534   char *cs;
4535   char **result;
4536
4537   result=vcalloc (len+1, sizeof (char*));
4538   array=sort_string_array (array, len);
4539
4540   for (c=0, a=0, n=0, cs=NULL; a< len; a++)
4541     {
4542       if (cs == NULL || !strm (cs, array[a]))
4543         {
4544           if (cs)
4545             {
4546               result[c]=vcalloc (strlen (cs)+20, sizeof (char));
4547               sprintf ( result[c], "%s %d", cs, n);
4548               c++;
4549             }
4550
4551           n=1;
4552           cs=array[a];
4553         }
4554       else n++;
4555     }
4556
4557   result[c]=vcalloc (strlen (cs)+20, sizeof (char));
4558   sprintf ( result[c], "%s %d", cs, n);
4559
4560   return result;
4561 }
4562
4563 int   get_first_non_white_char (char *name)
4564 {
4565   if ( !name) return 0;
4566   else if ( !check_file_exists (name))return 0;
4567   else
4568     {
4569       FILE *fp;
4570       char c;
4571
4572       fp=vfopen (name, "r");
4573       while (isspace ((c=fgetc (fp)))&& c!=EOF);
4574       vfclose (fp);
4575       return c;
4576     }
4577 }
4578 int   count_n_char_x_in_file(char *name, char x)
4579       {
4580       FILE *fp;
4581       int n, c;
4582
4583       n=0;
4584       fp=vfopen(name, "r");
4585       while ( (c=fgetc(fp))!=EOF)n+=(c==x);
4586       vfclose (fp);
4587       return n;
4588       }
4589 int   count_n_char_in_file(char *name)
4590       {
4591       int  c, n;
4592       FILE *fp;
4593
4594       n=0;
4595       fp=vfopen(name, "r");
4596       while ( (c=fgetc(fp))!=EOF){n++;}
4597       vfclose (fp);
4598       return n;
4599       }
4600 int count_n_line_in_file ( char *name )
4601      {
4602       int  c, n;
4603       FILE *fp;
4604
4605       n=0;
4606       fp=vfopen(name, "r");
4607       while ( (c=fgetc(fp))!=EOF)n+=(c=='\n');
4608       vfclose (fp);
4609       return n;
4610       }
4611 int measure_longest_line_in_file ( char *name )
4612      {
4613       int  c;
4614       FILE *fp;
4615       int longest=0;
4616       int current=0;
4617
4618
4619       fp=vfopen(name, "r");
4620       while ( (c=fgetc(fp))!=EOF)
4621           {
4622               if ( c=='\n'){longest=MAX(longest, current+1);current=0;}
4623               else current++;
4624           }
4625       longest=MAX(longest, current+1);
4626
4627       vfclose (fp);
4628       return longest;
4629       }
4630 char *input_name ()
4631         {
4632         char *string;
4633         int a;
4634         char ch;
4635
4636         string= vcalloc ( 500, sizeof ( char));
4637
4638         a=0;
4639         while ( ( ch=getchar())!='\n' && a<500)
4640                         string[a++]=ch;
4641         string[a]='\0';
4642
4643         if ( string[0]=='\0')
4644                 {
4645                 vfree (string);
4646                 return NULL;
4647                 }
4648         else
4649                 return string;
4650         }
4651
4652 FILE * vtmpfile()
4653     {
4654         return tmpfile();
4655     }
4656 int string_variable_isset (char *name)
4657 {
4658   if (store_string_variable (name,NULL, ISSET)==NULL)return 0;
4659   else return 1;
4660 }
4661 char* set_string_variable (char *name, char * v)
4662 {
4663   return store_string_variable (name, v, SET);
4664 }
4665 char * get_string_variable (char *name)
4666 {
4667   return store_string_variable (name,NULL, GET);
4668 }
4669 char* unset_string_variable (char *name)
4670 {
4671   return store_string_variable (name,0, UNSET);
4672 }
4673
4674 char* store_string_variable (char *name, char* v, int mode)
4675 {
4676   static char **name_array, **val_array;
4677   static int n;
4678   int a;
4679
4680   if ( mode == SET)
4681     {
4682
4683       for (a=0; a<n; a++)
4684         {
4685           if ( strm (name,name_array[a]))
4686             {
4687               if (v)
4688                 {
4689                   val_array[a]=vrealloc (val_array[a], strlen (v)+1);
4690                   sprintf (val_array[a],"%s",v);
4691                 }
4692               else val_array[a]='\0';
4693               return v;
4694             }
4695         }
4696       if (!name_array)
4697         {
4698           name_array=vcalloc (1, sizeof (char*));
4699           val_array=vcalloc  (1, sizeof (char*));
4700         }
4701       else
4702         {
4703           name_array=vrealloc (name_array, (n+1)*sizeof (char*));
4704           val_array=vrealloc (val_array, (n+1)*sizeof (char*));
4705         }
4706       name_array[n]=vcalloc ( strlen (name)+1, sizeof (char));
4707       val_array[n]=vcalloc ( strlen (v)+1, sizeof (char));
4708
4709       sprintf ( name_array[n], "%s", name);
4710       sprintf ( val_array[n], "%s", v);
4711       n++;
4712
4713       return v;
4714
4715     }
4716   else if ( mode == ISSET)
4717     {
4718        for (a=0; a<n; a++)
4719          if ( strm (name_array[a], name))return (char *)1;
4720     }
4721   else if ( mode == UNSET)
4722     {
4723       for (a=0; a<n; a++)
4724         if ( strm (name_array[a], name))
4725           {
4726             name_array[a][0]='\0';
4727             val_array[a][0]='\0';
4728             return 0;
4729           }
4730       add_warning (stdout, "WARNING: Could not UNSET the value of %s. You must SET the value before it is used", name);
4731     }
4732   else if (mode==GET)
4733     {
4734       for (a=0; a<n; a++)
4735         {
4736           if ( strm (name_array[a], name))
4737             return val_array[a];
4738         }
4739     }
4740   return NULL;
4741 }
4742 int int_variable_isset (char *name)
4743 {
4744   return store_int_variable (name,0, ISSET);
4745 }
4746 int set_int_variable (char *name, int v)
4747 {
4748   return store_int_variable (name, v, SET);
4749 }
4750 int get_int_variable (char *name)
4751 {
4752   return store_int_variable (name, 0, GET);
4753 }
4754 int unset_int_variable (char *name)
4755 {
4756   return store_int_variable (name,0, UNSET);
4757 }
4758
4759 int store_int_variable (char *name, int v, int mode)
4760 {
4761   static char **name_array;
4762   static int *val_array;
4763   static int n;
4764   int a;
4765
4766   if ( mode == SET)
4767     {
4768       for (a=0; a<n; a++)
4769         {
4770           if ( strm (name,name_array[a]))
4771             {
4772               val_array[a]=v;
4773               return v;
4774             }
4775         }
4776       if (!name_array)
4777         {
4778           name_array=vcalloc (1, sizeof (char*));
4779           val_array=vcalloc  (1, sizeof (int));
4780         }
4781       else
4782         {
4783           name_array=vrealloc (name_array, (n+1)*sizeof (char*));
4784           val_array=vrealloc (val_array, (n+1)*sizeof (int));
4785         }
4786       name_array[n]=vcalloc ( strlen (name)+1, sizeof (char));
4787       sprintf ( name_array[n], "%s", name);
4788       val_array[n]=v;
4789       n++;
4790
4791       return v;
4792
4793     }
4794   else if ( mode == ISSET)
4795     {
4796        for (a=0; a<n; a++)
4797          if ( strm (name_array[a], name))return 1;
4798     }
4799   else if ( mode == UNSET)
4800     {
4801       for (a=0; a<n; a++)
4802         if ( strm (name_array[a], name))
4803           {
4804             name_array[a][0]='\0';
4805             val_array[a]=0;
4806             return 0;
4807           }
4808       add_warning (stdout, "WARNING: Could not UNSET the value of %s. You must SET the value before it is used", name);
4809     }
4810   else if (mode==GET)
4811     {
4812       for (a=0; a<n; a++)
4813         {
4814           if ( strm (name_array[a], name))
4815             return val_array[a];
4816         }
4817     }
4818   return 0;
4819 }
4820
4821 char *add2t_coffee_environement (char *env, ...)
4822 {
4823   char *s;
4824   va_list ap;
4825   char *file;
4826   char *tmp_env_file;
4827   char *env_file;
4828   FILE *in, *out, *fp;
4829   static char *buf,*value1, *value2;
4830
4831   if (!buf)
4832     {
4833       buf=vcalloc (LONG_STRING, sizeof (char));
4834       value1=vcalloc (LONG_STRING, sizeof (char));
4835       value2=vcalloc (LONG_STRING, sizeof (char));
4836     }
4837   file=vtmpnam (NULL);
4838   tmp_env_file=vtmpnam (NULL);
4839
4840   va_start (ap, env);
4841   fp=vfopen (file, "w");
4842   vfprintf (fp,env, ap);
4843   vfclose (fp);
4844   va_end (ap);
4845   s=file2string (file);
4846   sprintf (value1, "%s", s);
4847   substitute_char (value1, '=', '\0');
4848
4849   //add the new environement variable to the env file
4850   if (!get_string_variable ("t_coffee_env"))get_t_coffee_environement (NULL);
4851   env_file=get_string_variable ("t_coffee_env");
4852
4853   out=vfopen(tmp_env_file, "w");
4854   if ( check_file_exists (env_file))
4855     {
4856       in=vfopen (env_file, "r");
4857
4858       while ( (fgets (buf,LONG_STRING, in))!=NULL)
4859         {
4860           sprintf (value2, "%s", buf);
4861           substitute_char (value2, '=', '\0');
4862           if ( !strm (value1, value2))fprintf ( out, "%s",buf);
4863         }
4864       vfclose (in);
4865     }
4866
4867   fprintf ( out, "%s\n", s);
4868   vfclose (out);
4869
4870   if ( check_file_exists (tmp_env_file))printf_system ("cp %s %s", tmp_env_file, env_file);
4871   //  fp=vfopen (get_string_variable ("t_coffee_env"), "a");
4872   //fprintf ( fp, "%s\n",s);
4873   //vfclose (fp);
4874
4875   //add the variable to T-Coffee
4876   putenv (s);
4877   return NULL;
4878 }
4879 char *get_t_coffee_environement (char *file)
4880 {
4881   char *def_env;
4882
4883
4884   if ( file);
4885   else if ((file=getenv ("ENV_4_TCOFFEE")));
4886   else if (get_dir_4_tcoffee())
4887     {
4888       file=vcalloc ( 100, sizeof (char));
4889       sprintf (file, "%s/t_coffee_env", get_dir_4_tcoffee());
4890     }
4891   else
4892     {
4893       file=vcalloc ( 100, sizeof (char));
4894       sprintf (file, "%s/t_coffee_env", getenv ("HOME"));
4895     }
4896   set_string_variable ("t_coffee_env", file);
4897
4898   def_env=vtmpnam (NULL);
4899   printf_system ("env > %s", def_env);
4900   file_putenv (def_env);
4901   file_putenv (file);
4902   if (!getenv ("DIR_4_TCOFFEE")) cputenv("DIR_4_TCOFFEE=%s", get_dir_4_tcoffee());
4903   if (!getenv ("TMP_4_TCOFFEE")) cputenv("TMP_4_TCOFFEE=%s", get_tmp_4_tcoffee());
4904   if (!getenv ("CACHE_4_TCOFFEE")) cputenv("CACHE_4_TCOFFEE=%s", get_cache_4_tcoffee());
4905
4906   set_path_4_plugins (NULL);
4907   get_email_from_env();
4908   get_proxy_from_env();
4909 }
4910 char * set_path_4_plugins (char *plugins)
4911 {
4912   if (plugins)get_plugins_4_tcoffee (plugins);
4913
4914   //set various important variables
4915   if (!getenv ("MAFFT_BINARIES"))cputenv ("MAFFT_BINARIES=%s", get_plugins_4_tcoffee(NULL));
4916   if (!getenv ("PLUGINS_4_TCOFFEE")) cputenv("PLUGINS_4_TCOFFEE=%s", get_plugins_4_tcoffee(NULL));
4917
4918   if (!getenv ("MCOFFEE_4_TCOFFEE")) cputenv("MCOFFEE_4_TCOFFEE=%s", get_mcoffee_4_tcoffee());
4919   if (!getenv ("METHODS_4_TCOFFEE")) cputenv("METHODS_4_TCOFFEE=%s", get_methods_4_tcoffee());
4920   cputenv ("PATH=%s:%s", get_plugins_4_tcoffee(NULL), getenv ("PATH"));
4921   //set various packages
4922   add_package2_tcoffee_env ("CLUSTALW_4_TCOFFEE");
4923   add_package2_tcoffee_env ("CLUSTALW2_4_TCOFFEE");
4924   add_package2_tcoffee_env ("FUGUE_4_TCOFFEE");
4925   add_package2_tcoffee_env ("CLUSTALW2_4_TCOFFEE");
4926   add_package2_tcoffee_env ("TMALIGN_4_TCOFFEE");
4927   add_package2_tcoffee_env ("SAP_4_TCOFFEE");
4928   add_package2_tcoffee_env ("DALILITEc_4_TCOFFEE");
4929   add_package2_tcoffee_env ("MUSTANG_4_TCOFFEE");
4930   add_package2_tcoffee_env ("NCBIBLAST_4_TCOFFEE");
4931   add_package2_tcoffee_env ("MAFFT_4_TCOFFEE");
4932   add_package2_tcoffee_env ("DIALIGNTX_4_TCOFFEE");
4933   add_package2_tcoffee_env ("POA_4_TCOFFEE");
4934   add_package2_tcoffee_env ("PROBCONS_4_TCOFFEE");
4935   add_package2_tcoffee_env ("PROBCONSRNA_4_TCOFFEE");
4936   add_package2_tcoffee_env ("TCOFFEE_4_TCOFFEE");
4937   add_package2_tcoffee_env ("PCMA_4_TCOFFEE");
4938   add_package2_tcoffee_env ("KALIGN_4_TCOFFEE");
4939   add_package2_tcoffee_env ("AMAP_4_TCOFFEE");
4940   add_package2_tcoffee_env ("PRODA_4_TCOFFEE");
4941   add_package2_tcoffee_env ("PRANK_4_TCOFFEE");
4942   add_package2_tcoffee_env ("CONSAN_4_TCOFFEE");
4943   add_package2_tcoffee_env ("RNAPlfold_4_TCOFFEE");
4944   add_package2_tcoffee_env ("HMMtop_4_TCOFFEE");
4945   add_package2_tcoffee_env ("GOR4_4_TCOFFEE");
4946   return NULL;
4947 }
4948
4949 int add_package2_tcoffee_env (char *package)
4950 {
4951   char *v1, *v2;
4952
4953   if ((v1=getenv (package))==NULL)return 0;
4954
4955   v2=filename2path (v1);
4956   if (!v2){vfree (v1);return 0;}
4957   cputenv ("PATH=%s:%s",v2, getenv ("PATH"));
4958   vfree (v1); vfree (v2);
4959   return 1;
4960 }
4961
4962
4963 char * Proxy( int input_mode, int write_mode)
4964 {
4965   static char *proxy;
4966   static int set;
4967
4968   if ( write_mode==RESET){set=0;vfree (proxy); proxy=NULL;return Proxy(input_mode, SET);}
4969
4970   if (set);
4971   else if ( input_mode ==INPUT){proxy=get_proxy();}
4972   else if ( input_mode ==ENV  ){proxy=get_proxy_from_env (); }
4973   else printf_exit (EXIT_FAILURE,stderr, "Unknown mode for Proxy [FATAL:%s(%s)]", PROGRAM, VERSION);
4974
4975   set=1;
4976   if (proxy)set_proxy (proxy);
4977   else proxy=vcalloc (1, sizeof (char));
4978
4979   return proxy;
4980 }
4981 char * get_proxy_from_env ()
4982 {
4983   char *proxy=NULL;
4984
4985   if ((proxy=get_string_variable ("cl_proxy"))){;}//Command line proxy always wins
4986   else if ((proxy=getenv ("http_proxy_4_TCOFFEE")));
4987   else if ((proxy=get_string_variable ("proxy")));//use default T-Coffee proxy
4988   else if ( getenv ("HTTP_proxy") && getenv ("http_proxy")){return getenv ("HTTP_proxy");}//get environement proxy
4989   else if ((proxy=getenv ("HTTP_proxy")));//id
4990   else if ((proxy=getenv ("http_proxy")));//id
4991   else if ((proxy=getenv ("HTTP_PROXY")));//id
4992   else if ((proxy=getenv ("ALL_proxy")));//id
4993   else if ((proxy=getenv ("all_proxy")));//id
4994   else if ((proxy=getenv ("ALL_PROXY")));//id
4995
4996   if (proxy)set_proxy(proxy);
4997   return proxy;
4998 }
4999 char *get_proxy ()
5000 {
5001
5002   char *proxy=NULL;
5003
5004   if ( (proxy=get_proxy_from_env()) && simple_check_internet_connection(NULL));
5005   else //read in the proxy
5006     {
5007       proxy_message();
5008       input_proxy();
5009       proxy=get_proxy();
5010     }
5011   return proxy;
5012 }
5013 int set_proxy (char *proxy)
5014 {
5015   char *http;
5016   if (!proxy) return 0;
5017
5018   cputenv ("HTTP_proxy_4_TCOFFEE=%s", proxy);
5019   cputenv ("HTTP_proxy=%s", proxy);
5020   cputenv ("http_proxy=%s", proxy);
5021   cputenv ("HTTP_PROXY=%s", proxy);
5022   cputenv ("ALL_proxy=%s", proxy);
5023   cputenv ("ALL_PROXY=%s", proxy);
5024   cputenv ("all_proxy=%s", proxy);
5025
5026   return 1;
5027 }
5028 char *proxy_message ()
5029 {
5030   fprintf ( stderr, "\n\n");
5031   fprintf ( stderr, "*************************************************************************************************\n");
5032   fprintf ( stderr, "*                        IMPORTANT: Please Read Carefuly                                        *\n");
5033   fprintf ( stderr, "*                                                                                               *\n");
5034   fprintf ( stderr, "* ------------------------------FIRST RUN CONFIGURATION---------------------------------------- *\n");
5035   fprintf ( stderr, "*                                                                                               *\n");
5036   fprintf ( stderr, "*                                                                                               *\n");
5037   fprintf ( stderr, "* If you are behind a firewall, you must enter your proxy address to use webservices            *\n");
5038   fprintf ( stderr, "* This address is usualy something like: http://some.place.here:8080                            *\n");
5039   fprintf ( stderr, "*                                                                                               *\n");
5040   fprintf ( stderr, "* The proxy you will provide can be changed anytime by editing the file:                        *\n");
5041   fprintf ( stderr, "*  %s/t_coffee_env                                   \n",get_dir_4_tcoffee());
5042   fprintf ( stderr, "*************************************************************************************************\n");
5043   return NULL;
5044 }
5045 char *input_proxy ()
5046 {
5047   char *proxy;
5048   char *answer;
5049   static int ntries;
5050
5051   ntries++;
5052   if ( ntries==MAX_N_TRIES)
5053     printf_exit (EXIT_FAILURE,stderr, "\nERROR: Could not use provided proxy. Provide proxy via the -email flag [FATAL:%s]\n", PROGRAM);
5054
5055   fprintf ( stderr, "\n ### Tip: If you work from home (ADSL) you probably do no need a proxy. ");
5056   fprintf ( stderr, "\n ### Tip: The proxy is the address you may have had to enter in your navigator. ");
5057
5058   fprintf ( stderr, "\n ### Proxy Input, Enter your proxy and type Return (Type Return if you do not need a proxy): ");
5059   proxy=input_name();
5060
5061   if ( !proxy)
5062     {
5063       proxy=vcalloc ( 100, sizeof (char));
5064       sprintf ( proxy, "");
5065     }
5066   else if (!strstr (proxy, "."))return input_proxy();
5067   else
5068     {
5069       int a, c;
5070       a=0;
5071       while ((c=proxy[a])!='\0')
5072         {
5073           if (!isalnum(c) && c!='.' && c!='_' )return input_proxy();
5074           a++;
5075         }
5076     }
5077   fprintf ( stderr, "\n\tYou have entered:[%s]%s \nIs this correct ([Y] or N)?: ",proxy, (strm (proxy, "")?"(No Proxy Needed)":""));
5078   answer=input_name();
5079
5080   if ( !answer || answer[0]=='y' ||answer[0]=='Y')
5081     {
5082       vfree (answer);
5083     }
5084   else
5085     {
5086       vfree (answer);
5087       return input_proxy ();
5088     }
5089   add2t_coffee_environement ("http_proxy_4_TCOFFEE=%s", proxy);
5090   set_string_variable ("proxy", proxy);
5091   return proxy;
5092 }
5093
5094
5095
5096 char *input_email ()
5097 {
5098   char *email;
5099   char *answer;
5100   static int ntries;
5101
5102   ntries++;
5103
5104   if ( ntries==MAX_N_TRIES)
5105     printf_exit (EXIT_FAILURE, stderr,"\nERROR: Could not use provided Email. Provide Email via the -email flag [FATAL:%s]\n", PROGRAM);
5106
5107
5108
5109   fprintf ( stderr, "\n ### EMAIL Input, Enter your email and type Return: ");
5110   email=input_name();
5111   fprintf ( stderr, "\n\tYou have entered:%s \nIs this correct ([Y] or N)?: ", email);
5112   answer=input_name();
5113
5114   if ( !answer || answer[0]=='y' ||answer[0]=='Y')
5115     {
5116       vfree (answer);
5117     }
5118   else
5119     {
5120       vfree (answer);
5121       return input_email ();
5122     }
5123   add2t_coffee_environement("EMAIL_4_TCOFFEE=%s", email);
5124   set_string_variable ("email", email);
5125   return email;
5126 }
5127 static char *ebi_email_message ();
5128 char *ebi_email_message ()
5129 {
5130   fprintf ( stderr, "\n\n");
5131   fprintf ( stderr, "*************************************************************************************************\n");
5132   fprintf ( stderr, "*                        IMPORTANT: Please Read Carefuly                                        *\n");
5133   fprintf ( stderr, "*                                                                                               *\n");
5134   fprintf ( stderr, "* ------------------------------FIRST RUN CONFIGURATION---------------------------------------- *\n");
5135   fprintf ( stderr, "*                                                                                               *\n");
5136   fprintf ( stderr, "* Some commands of T-Coffee use the  EBI BLAST webservices. The EBI requires a valid E-mail     *\n");
5137   fprintf ( stderr, "* address for this service to be used (check: www.ebi.ac.uk/Tools/webservices/).                *\n");
5138   fprintf ( stderr, "* T-Coffee will keep it for further run but only to be used with the EBI webservices.           *\n");
5139   fprintf ( stderr, "*                                                                                               *\n");
5140   fprintf ( stderr, "* !!!!!!!!!!!!!!!!!! Your Email will not be sent to us, ONLY to the EBI !!!!!!!!!!!!!!!!!!!!!!!!*\n");
5141   fprintf ( stderr, "*                                                                                               *\n");
5142   fprintf ( stderr, "* -blast_server=EBI is the default mode of T-Coffee. If you do NOT want to provide your E-mail  *\n");
5143   fprintf ( stderr, "* you can use:                                                                                  *\n");
5144   fprintf ( stderr, "*    -blast_server=NCBI     (NCBI netblast)                                                     *\n");
5145   fprintf ( stderr, "*    -blast_server=LOCAL     Local NCBI BLAST                                                   *\n");
5146   fprintf ( stderr, "*                                                                                               *\n");
5147   fprintf ( stderr, "*The address you provide can be changed anytime by editing the file:                            *\n");
5148   fprintf ( stderr, "*    %s/t_coffee_env                                   \n",get_dir_4_tcoffee());
5149   fprintf ( stderr, "*************************************************************************************************\n");
5150   return NULL;
5151 }
5152 char *Email( int input_mode, int write_mode)
5153 {
5154   static char *email;
5155   static int set;
5156
5157   if ( write_mode==RESET){set=0;vfree (email); email=NULL;return Email(input_mode, SET);}
5158
5159   if (set);
5160   else if (input_mode ==INPUT){email=get_email();}
5161   else if (input_mode ==ENV  ){email=get_email_from_env (); }
5162   else printf_exit (EXIT_FAILURE, stderr,"Unknown mode for Email [FATAL:%s]",PROGRAM);
5163
5164   set=1;
5165
5166   if ( email)set_email (email);
5167   else email=vcalloc (1, sizeof (char));
5168   return email;
5169 }
5170 char *get_email_from_env ()
5171 {
5172   char *email;
5173
5174   if ( (email=get_string_variable ("cl_email")));
5175   else if ( (email=get_string_variable ("email")));
5176   else if ( (email=getenv ("EMAIL_4_TCOFFEE")));
5177   else if ( (email=getenv ("EMAIL")));
5178   return email;
5179 }
5180 char *get_email ()
5181 {
5182   char file [1000];
5183   static char *email;
5184
5185   if ( (email=get_email_from_env ()));
5186   else
5187     {
5188       ebi_email_message();
5189       email=input_email();
5190     }
5191
5192   if (!strstr (email, "@"))
5193     {
5194       ebi_email_message();
5195       email=input_email();
5196     }
5197   set_email (email);
5198
5199   return email;
5200 }
5201 int set_email (char *email)
5202 {
5203   if (!email) return 0;
5204
5205   cputenv ("EMAIL_4_TCOFFEE=%s", email);
5206   cputenv ("EMAIL=%s",email);
5207   return 1;
5208 }
5209 char *chomp (char *name)
5210 {
5211   int a=0;
5212   while ( name[a]!='\n' && name[a]!='\0')a++;
5213   name[a]='\0';
5214   return name;
5215 }
5216 static Tmpname *tmpname;
5217 static Tmpname *ntmpname;
5218
5219 static int n_tmpname;
5220 static int file2remove_flag;
5221
5222 char *set_file2remove_extension (char *extension, int mode)
5223 {
5224   static char ext[100];
5225   if (mode==SET)sprintf (ext, "%s", extension);
5226   else if ( mode==UNSET) ext[0]='\0';
5227   else if ( mode==GET);
5228   return ext;
5229 }
5230 int flag_file2remove_is_on ()
5231 {
5232   return file2remove_flag;
5233 }
5234 void set_file2remove_on()
5235 {
5236   file2remove_flag=1;
5237 }
5238 void set_file2remove_off()
5239 {
5240   file2remove_flag=0;
5241 }
5242
5243 char *add2file2remove_list (char *name)
5244 {
5245
5246
5247   if ( !tmpname || !name)ntmpname=tmpname=vcalloc ( 1, sizeof (Tmpname));
5248   else if (!ntmpname->name);
5249   else ntmpname=ntmpname->next=vcalloc ( 1, sizeof (Tmpname));
5250
5251   if (!name) return NULL;
5252
5253   ntmpname->name=vcalloc(strlen(name)+1, sizeof (char));
5254
5255   sprintf (ntmpname->name, "%s", name);
5256   return ntmpname->name;
5257 }
5258 //char *short_tmpnam_2(char *s);//used to generate very compact tmp names
5259 void  initiate_vtmpnam (char *file)
5260 {
5261   add2file2remove_list (NULL);
5262   tmpnam_2(NULL);
5263 }
5264 char *vtmpnam ( char *s1)
5265 {
5266   char *s,*s2;
5267
5268   n_tmpname++;
5269
5270   standard_initialisation(NULL, NULL);
5271
5272   s=vcalloc ( VERY_LONG_STRING, sizeof (char));
5273   s[0]='\0';
5274
5275   s=tmpnam_2 (s);
5276
5277   s2=add2file2remove_list (s);
5278   if (s!=s2)vfree (s);
5279   if (s1){sprintf (s1, "%s",s2);return s1;}
5280   else return s2;
5281 }
5282
5283
5284
5285 int get_vtmpnam2_root()
5286 {
5287   int MAX_TMPNAM_ROOT=10000;
5288   static int v;
5289
5290   if (v) ;
5291   else
5292     {
5293       vsrand(0);
5294       v=rand()%MAX_TMPNAM_ROOT;
5295     }
5296   return v;
5297 }
5298 char *tmpnam_2 (char *s)
5299      {
5300        static int root;
5301        static int file;
5302        char buf[VERY_LONG_STRING];
5303        static char root2[VERY_LONG_STRING];
5304        static char *tmpdir;
5305        static int name_size;
5306
5307        if ( !root || !s)
5308          {
5309            char *vtmpnam_prefixe;
5310
5311            name_size=MAX( 2*L_tmpnam, MAXNAMES*2)+1;
5312            root=get_vtmpnam2_root();
5313            sprintf ( root2, "%d%d_", root, (int)getpid());
5314
5315            vtmpnam_prefixe=vcalloc (strlen (root2)+strlen (get_tmp_4_tcoffee())+2, sizeof (char));
5316            sprintf (vtmpnam_prefixe, "%s/%s", get_tmp_4_tcoffee(), root2);
5317            set_string_variable ("vtmpnam_prefixe1", vtmpnam_prefixe);
5318            set_string_variable ("vtmpnam_prefixe2", root2);
5319            vfree (vtmpnam_prefixe);
5320          }
5321
5322        if (!s)return NULL;
5323        tmpdir=get_tmp_4_tcoffee();
5324
5325        sprintf (buf, "%s/%s%d_TCtmp%s",tmpdir,root2, file++,set_file2remove_extension (NULL, GET));
5326        if ( strlen(buf)>=name_size)s=vrealloc (s,(strlen(buf)+1)*sizeof (char));
5327        sprintf (s, "%s", buf);
5328        return s;
5329      }
5330 char *short_tmpnam_2(char *s)
5331 {
5332   static int root;
5333   static int file;
5334   char buf[VERY_LONG_STRING];
5335   static char root2[VERY_LONG_STRING];
5336   static char *tmpdir;
5337   static int name_size;
5338
5339   if ( !root || !s)
5340     {
5341       char *vtmpnam_prefixe;
5342
5343       name_size=MAX( 2*L_tmpnam, MAXNAMES*2)+1;
5344       root=get_vtmpnam2_root();
5345       sprintf ( root2, "%d%d", root,getpid());
5346
5347       vtmpnam_prefixe=vcalloc (strlen (root2)+strlen (get_tmp_4_tcoffee())+2, sizeof (char));
5348       sprintf (vtmpnam_prefixe, "%s", root2);
5349       set_string_variable ("vtmpnam_prefixe1", vtmpnam_prefixe);
5350       set_string_variable ("vtmpnam_prefixe2", root2);
5351       vfree (vtmpnam_prefixe);
5352     }
5353   if (!s) return NULL;
5354
5355   sprintf (buf, "%s%d%s",root2, file++,set_file2remove_extension (NULL, GET));
5356   if ( strlen(buf)>=name_size)s=vrealloc (s,(strlen(buf)+1)*sizeof (char));
5357   sprintf (s, "%s", buf);
5358   return s;
5359 }
5360
5361 char *vremove2 (char *s)
5362 {
5363   char command[1000];
5364   char list_file[1000];
5365   char ***list;
5366   int a;
5367
5368
5369   //Remove filenames with a wildcard
5370
5371   sprintf (list_file, "list_file_%d", (int)getpid());
5372   sprintf (command, "ls -1 %s>%s 2>/dev/null", s, list_file);
5373   safe_system (command);
5374
5375   list=file2list (list_file, " ");
5376
5377   a=0;
5378   while (list && list[a])
5379     {
5380       if ( check_file_exists (list[a][1]))
5381         {
5382           vremove (list[a][1]);
5383         }
5384       a++;
5385     }
5386   vremove (list_file);
5387   return NULL;
5388 }
5389 char *vremove (char *s)
5390 {
5391
5392
5393   if ( s && strstr (s, "*"))return vremove2(s);
5394   else if ( !s || !check_file_exists(s) ) return NULL;
5395   else if ( isdir (s))
5396     {
5397       rmdir (s);
5398       return NULL;
5399     }
5400   else
5401     {
5402       remove (s);
5403       return NULL;
5404     }
5405   return NULL;
5406 }
5407 int log_function ( char *fname)
5408 {
5409   char command[1000];
5410
5411   if ( check_file_exists (error_file))
5412     {
5413
5414       sprintf ( command, "cp %s %s", error_file, fname);
5415       my_system ( command);
5416       fprintf( stderr,"\n\t******************************************************************");
5417       fprintf( stderr, "\n\t* Full Log of [%s, %s] in File [%s]",PROGRAM, VERSION, fname);
5418       fprintf( stderr, "\n\t******************************************************************\n");
5419     }
5420   return 1;
5421 }
5422
5423 void clean_exit ()
5424 {
5425
5426   myexit (global_exit_signal);
5427 }
5428 void error_exit ( )
5429      {
5430        char command[1000];
5431        char final_report[1000];
5432
5433        if ( no_error_report)return;
5434        if ( check_file_exists (error_file))
5435          {
5436            sprintf ( final_report, "error_report.%s",PROGRAM);
5437            if ( !getenv ("NO_ERROR_REPORT_4_TCOFFEE"))sprintf ( command, "cp %s %s", error_file, final_report);
5438            my_system ( command);
5439            fprintf( stderr,"\n\t******************************************************************");
5440            fprintf( stderr, "\n\t* Job NOT Completed:[%s, %s]",PROGRAM, VERSION);
5441            fprintf( stderr, "\n\t* Please CHECK:                                       ");
5442            fprintf( stderr, "\n\t* \t-1 The format of your Input Files                 ");
5443            fprintf( stderr, "\n\t* \t-2 The parameters                                 ");
5444            fprintf( stderr, "\n\t* \t-3 The use of special characters in sequence names:");
5445            fprintf( stderr, "\n\t* \t\t (@, |, %%...)");
5446
5447            fprintf( stderr, "\n\t* \t-4 The Online Doc (%s)                   ", URL);
5448            fprintf( stderr, "\n\t* \t-5 Send the file %s to:", final_report);
5449            fprintf( stderr, "\n\t* \t\t%s",EMAIL);
5450
5451            fprintf( stderr, "\n\t* If you run T-Coffee over the WEB:");
5452            fprintf( stderr, "\n\t* \tWindows Cut and Paste is sometimes erratic and");
5453            fprintf( stderr, "\n\t* \tit can loose carriage returns. If you suspect this,");
5454            fprintf( stderr, "\n\t* \ttry to cut and paste through an intermediate application");
5455            fprintf( stderr, "\n\t* \t(word pad) and inspect the results\n\n");
5456            fprintf( stderr, "\n\t* CONFIDENTIALITY:");
5457            fprintf( stderr, "\n\t* \tThe File %s may contain your personnal DATA", final_report);
5458            fprintf( stderr, "\n\t* \tRemove ALL confidential DATA from this file BEFORE sending it");
5459            fprintf( stderr, "\n\t******************************************************************\n");
5460            print_command_line(stderr);
5461          }
5462        print_exit_failure_message ();
5463        myexit (global_exit_signal=EXIT_FAILURE);
5464      }
5465 void main_exit ()
5466 {
5467
5468   clean_function ();
5469   exit (global_exit_signal);
5470 }
5471
5472
5473 void clean_function ()
5474 {
5475          Tmpname *b;
5476          char *tmp;
5477          static int done;
5478          Tmpname *start;
5479
5480          start=tmpname;
5481
5482          if ( done==1) return;
5483          else done=1;
5484
5485          kill_child_pid();
5486          add_method_output2method_log (NULL, NULL, NULL, NULL, decode_name (NULL, CODELIST));
5487
5488          //Kill all child processes (if any)
5489
5490
5491          if (getenv ("DEBUG_TMP_FILE") && atoi (getenv("DEBUG_TMP_FILE"))==1)
5492            {
5493              fprintf ( stderr, "\n[DEBUG_TMP_FILE:%s] TEMPORARY FILES HAVE NOT Been Removed:", PROGRAM);
5494              while (start)
5495                  {
5496                    if ( getenv("PRINT_TMPFILE_NAME") && atoi(getenv("PRINT_TMPFILE_NAME"))==1)
5497                      {
5498                        fprintf ( stderr, "\n\t%s [EXISTS:%s]", tmpname->name, (check_file_exists(tmpname->name))?"YES":"NO");
5499                      }
5500
5501                    b=start;
5502                    start=start->next;
5503                    vfree(b->name);
5504                    vfree(b);
5505                  }
5506              fprintf ( stderr, "\n");
5507              return;
5508            }
5509          else
5510            {
5511              char name[10000];
5512              char *x;
5513
5514              //printf_system ("clean_cache.pl -dir=%s -size=0 -age=10 -force",get_tmp_4_tcoffee ());
5515              //Remove all the temporary files generated during THIS run
5516              while (start)
5517                {
5518
5519                  remove (start->name);
5520                  if (isdir(start->name))rrmdir (start->name);
5521
5522                  b=start;
5523                  start=start->next;
5524                  vfree(b->name);vfree(b);
5525                  }
5526
5527              //Remove potential log
5528              //remove ( TO_NULL_DEVICE);
5529              //Remove all the tmp/tmp* files generated by this process [potentialy generated by secondary processes]
5530              //if ((x=get_string_variable("vtmpnam_prefixe1"))){ sprintf ( name, "%s*", x); vremove(name);}
5531              //Remove all the tmp* file that may have sneaked in (from clustalw for instance) and could be in the current dir
5532              //if ((x=get_string_variable("vtmpnam_prefixe2"))){ sprintf ( name, "%s*", x); vremove(name);}
5533              //remove (get_tmp_4_tcoffee ());//removes tmpdir if empty
5534            }
5535          return;
5536      }
5537 FILE *NFP;/*Null file pointer: should only be open once*/
5538
5539 /*********************************************************************/
5540 /*                                                                   */
5541 /*                         CACHE_FUNCTION                            */
5542 /*                                                                   */
5543 /*                                                                   */
5544 /*********************************************************************/
5545 static char *cache;
5546 char * prepare_cache ( const char *mode)
5547 {
5548   char command[1000];
5549
5550   cache =vcalloc ( 10000, sizeof(char));
5551
5552   if (strm (mode, "use"))
5553     {
5554       sprintf (cache, "%s",get_cache_4_tcoffee());
5555     }
5556
5557   else if ( strm (mode, "ignore") ||  strm (mode, "no"))
5558     {
5559
5560       cache=vtmpnam(cache);
5561       strcat (cache, "/");
5562       sprintf ( command, "mkdir %s",cache);
5563       my_system ( command);
5564     }
5565   else if ( strm (mode, "update"))
5566     {
5567       cache=vtmpnam(cache);
5568       strcat (cache, "/");
5569       sprintf ( command, "mkdir %s",cache);
5570       my_system ( command);
5571     }
5572   else if ( strm (mode, "local"))
5573     {
5574       cache[0]='\0';
5575     }
5576   else
5577     {
5578       sprintf ( cache, "%s/",mode);
5579       my_mkdir ( cache);
5580     }
5581   return cache;
5582
5583 }
5584
5585 char * get_cache_dir()
5586 {
5587   if ( cache==NULL){cache=vcalloc (1, sizeof (char));cache[0]='\0';}
5588   return cache;
5589 }
5590
5591 void update_cache ()
5592 {
5593   char command[1000];
5594   char old_cache[1000];
5595
5596   sprintf ( old_cache, "%s", get_cache_dir());
5597   prepare_cache( "use");
5598   sprintf ( command, "mv %s* %s",old_cache, get_cache_dir());
5599   my_system (command);
5600   sprintf ( command, "rmdir %s",old_cache);
5601   my_system (command);
5602 }
5603 void ignore_cache()
5604 {
5605   char command[1000];
5606
5607   if (getenv4debug ("DEBUG_TMP_FILE"))
5608     {
5609       fprintf ( stderr, "\n[DEBUG_TMP_FILE:%s] TEMPORARY CACHE HAS NOT Been Removed:\n\t%s\n", PROGRAM,get_cache_dir());
5610     }
5611   else
5612     {
5613
5614       sprintf ( command, "rm -r %s",get_cache_dir());
5615       my_system (command);
5616     }
5617   return;
5618
5619 }
5620
5621
5622 FILE * vfopen  ( char *name_in, char *mode)
5623     {
5624     FILE *fp;
5625     int get_new_name;
5626     int tolerate_mistake;
5627     int cache_used=0;
5628     FILE *tmp_fp;
5629     int c;
5630     static char *name;
5631     static char *name2;
5632     static char *stdin_file;
5633
5634
5635     if ( !name_in)return NULL;
5636     if (!name){name=vcalloc (1000, sizeof (char));}
5637     if (!name2){name2=vcalloc (1000, sizeof (char));}
5638
5639     sprintf ( name, "%s", name_in);
5640     tild_substitute (name, "~", get_home_4_tcoffee());
5641
5642     get_new_name=tolerate_mistake=0;
5643     if ( mode[0]=='g'){get_new_name=1; mode++;}
5644     else if ( mode[0]=='t'){tolerate_mistake=1;mode++;}
5645 /*Use the cached version from CACHE_4_TCOFFEE*/
5646     else if ( mode[0]=='c'){cache_used=1;mode++;}
5647
5648     if (name==NULL ||strm5 ( name, "no","NO","No","NULL","/dev/null") || strm2 (name, "no_file", "NO_FILE"))
5649                 {
5650                 if ( NFP==NULL)NFP=fopen (NULL_DEVICE, mode);
5651                 return NFP;
5652                 }
5653     else if ( strm3 (name,"stderr","STDERR","Stderr"))return stderr;
5654     else if ( strm3 (name,"stdout","STDOUT","Stdout"))return stdout;
5655     else if ( strm3 ( name, "stdin","STDIN","Stdin"))
5656         {
5657           if (!stdin_file)
5658             {
5659               stdin_file=vtmpnam (NULL);
5660               tmp_fp=vfopen ( stdin_file, "w");
5661               while ( (c=fgetc(stdin))!=EOF)fprintf (tmp_fp, "%c", c);
5662               vfclose ( tmp_fp);
5663             }
5664           return vfopen (stdin_file, "r");
5665         }
5666
5667     else if ( strm (name, "") && (strm (mode, "w") ||strm (mode, "a")) )return stdout;
5668     else if ( strm (name, "") && strm (mode, "r"))return stdin;
5669     else if ( (fp= fopen ( name, mode))==NULL)
5670                 {
5671                 if ( strcmp (mode, "r")==0 && cache_used==0)
5672                   {
5673                     sprintf ( name2, "%s%s",get_cache_dir(), name);
5674                     return vfopen ( name2, "cr");
5675                   }
5676                 else if ( strcmp (mode, "r")==0 && cache_used==1)
5677                         {
5678                         fprintf (stderr, "\nCOULD NOT READ %s\n", name);
5679                         if ( get_new_name){fprintf ( stderr, "\nNew name: ");return vfopen (input_name(), mode-1);}
5680                         else if ( tolerate_mistake)return NULL;
5681                         else
5682                             {
5683                             fprintf (stderr, "\nFORCED EXIT (NON INTERACTIVE MODE)\n");
5684                             if ( getenv ( "DEBUG_TCOFFEE"))crash ("DEBUG");
5685                             else myexit (EXIT_FAILURE);
5686                             }
5687                         }
5688                 else if ( strcmp (mode, "a")==0 && cache_used==0)
5689                   {
5690                     sprintf ( name2, "%s%s",get_cache_dir(), name);
5691                     return vfopen ( name, "ca");
5692                   }
5693                 else if ( strcmp (mode, "a")==0 && cache_used==1)
5694                         {
5695                         fprintf (stderr, "\nCOULD NOT Append anything to %s\n", name);
5696                         if ( get_new_name){fprintf ( stderr, "\nNew name: ");return vfopen (input_name(), mode-1);}
5697                         else if ( tolerate_mistake)return NULL;
5698                         else
5699                             {
5700                             fprintf (stderr, "\nFORCED EXIT (NON INTERACTIVE MODE)\n");
5701                             myexit (EXIT_FAILURE);
5702                             }
5703                         }
5704                 else if ( strcmp (mode, "w")==0)
5705                         {
5706                         fprintf (stderr, "\nCANNOT WRITE %s\n", name);
5707                         if ( get_new_name==1){fprintf ( stderr, "\nNew name: ");return vfopen (input_name(), mode-1);}
5708                         else if ( tolerate_mistake)return NULL;
5709                         else
5710                             {
5711                             fprintf (stderr, "\nFORCED EXIT (NON INTERACTIVE MODE): %s %s\n", (strcmp ( mode, "r")==0)?"READ":"WRITE", name);
5712                             myexit(EXIT_FAILURE);
5713                             }
5714                         }
5715                 }
5716     else
5717         return fp;
5718
5719     return NULL;
5720     }
5721
5722 FILE * vfclose ( FILE *fp)
5723        {
5724        if ( fp==NFP)return NULL;
5725        if ( fp==stdout)return stdout;
5726        if ( fp==stderr)return stderr;
5727        if ( fp==stdin) return stdin;
5728        if ( fp==NULL)return NULL;
5729        else fclose (fp);
5730        return NULL;
5731        }
5732
5733
5734 int echo ( char *string, char *fname)
5735 {
5736 int a;
5737 /*
5738 description:
5739 prints the content of string into file fname
5740
5741 in:
5742 string= string to print
5743 fname =name of the file to create
5744 */
5745
5746 FILE *fp;
5747
5748     fp=vfopen ( fname, "w");
5749     fprintf (fp, "%s", string);
5750     a=fclose (fp);
5751     return a;
5752
5753 }
5754
5755 int   file_cat ( char *from, char *to)
5756 {
5757   FILE *fp;
5758   //appends the content of file1 to file 2
5759   if (!(fp=vfopen (to, "a")))return 0;
5760   if (!display_file_content (fp, from)) return 0;
5761   vfclose (fp);
5762   return 1;
5763 }
5764
5765 FILE* display_file_content (FILE *output, char *name)
5766 {
5767   FILE *fp;
5768   int c;
5769   if ( !name || !check_file_exists (name) || !(fp=vfopen (name, "r")))return NULL;
5770   while ( (c=fgetc(fp))!=EOF)fprintf (output,"%c", c);
5771   vfclose (fp);
5772   return output;
5773 }
5774
5775 char ***file2list ( char *name, char *sep)
5776 {
5777   /*Rturns an array where
5778     list[0]: first line
5779     list[0][0]: number of words
5780     list[0][1]:first word;
5781     list[n]=NULL
5782   */
5783   char **lines, ***list;
5784   int a, n;
5785
5786   lines=file2lines (name);
5787   if (!lines) return NULL;
5788   else
5789     {
5790       n=atoi (lines[0]);
5791
5792       list=vcalloc ( n+1, sizeof (char**));
5793       for ( a=1; a<n; a++)
5794         {
5795
5796           list[a-1]=string2list2(lines[a], sep);
5797         }
5798     }
5799
5800   free_arrayN((void**)lines, 2);
5801   return list;
5802 }
5803 char **file2lines (char *name)
5804 {
5805   /*lines[0]->nlines;
5806     lines[1]->first_line
5807   */
5808   char **lines;
5809   char *string;
5810
5811
5812   string=file2string (name);
5813   if ( !string) return NULL;
5814   else
5815     {
5816       lines=string2list2(string, "\n");
5817       vfree ( string);
5818       return lines;
5819     }
5820 }
5821
5822 char *string2file ( char *string, char *file, char *mode)
5823 {
5824   FILE *fp;
5825   if (!file)
5826     file=vtmpnam (NULL);
5827   fp=vfopen (file, mode);
5828   fprintf (fp, "%s", string);
5829   vfclose (fp);
5830   return file;
5831 }
5832 char *file2string (char *name)
5833 {
5834   FILE*fp;
5835   char *string;
5836   int a, c;
5837
5838   if (!name || !check_file_exists (name))return NULL;
5839   else
5840     {
5841       string=vcalloc ( count_n_char_in_file(name)+1, sizeof (char));
5842       fp=vfopen (name, "r");a=0;
5843       while ( (c=fgetc(fp))!=EOF)
5844         {
5845           string[a++]=c;
5846         }
5847       string[a]='\0';
5848       vfclose (fp);
5849       return string;
5850     }
5851 }
5852
5853 int get_cl_param (int argc, char **argv, FILE **fp,char *para_name, int *set_flag, char *type, int optional, int max_n_val,char *usage, ...)
5854         {
5855         /*
5856         usage:
5857         argc:       n_ arg
5858         argv        list *
5859         para_name   param
5860         set_flag    set to 1 if param set;
5861         para_type   F, I, S, R_FN (read_file, name), W_FN (written file, name), R_FP (pointer)
5862         max_n_val   maximum number of values;
5863         optional    1 for yes, 0 for no
5864         usage       usage list with optional value;
5865         val         pointer to the varaible holding the value(s)
5866         default1    default value (if value id not there)
5867         default2    default value if the flag is there but no value set ("")indicates an error
5868         range_left  min value ( "any" for any);
5869         range_right  max_value ( "any" for any);
5870         */
5871         int pos=0;
5872         int a;
5873         va_list ap;
5874
5875         int   *int_val=NULL;
5876         float *float_val=NULL;
5877         char  **string_val=NULL;
5878
5879
5880         char  *range_right;
5881         char  *range_left;
5882
5883
5884         char *default_value1;
5885         char *default_value2;
5886         int n_para=0;
5887         double max, min;
5888
5889         static char **parameter_list;
5890         static int    number_of_parameters;
5891
5892         char **para_name_list;
5893         int    n_para_name;
5894
5895         char **para_val;
5896         int    n_para_val;
5897
5898         char **pv_l=NULL;
5899         int    n_pv_l;
5900         char **pv_r=NULL;
5901         int    n_pv_r;
5902         char   value[STRING];
5903
5904
5905
5906 /*CHECK THAT ALL THE PARAM IN ARG EXIST*/
5907         if ( para_name==NULL)
5908            {
5909            for ( a=1; a< argc; a++)
5910                {
5911                if ( is_parameter ( argv[a]))
5912                    if ( name_is_in_list ( argv[a], parameter_list, number_of_parameters, STRING)==-1)
5913                       {
5914                         fprintf ( stderr, "\n%s IS NOT A PARAMETER  OF %s [FATAL/%s %s]\n",argv[a], argv[0], argv[0], VERSION);
5915                       myexit(EXIT_FAILURE);
5916                       }
5917
5918
5919                }
5920
5921            free_char (parameter_list,-1);
5922            return 0;
5923            }
5924
5925         if ( parameter_list==NULL)parameter_list=declare_char(MAX_N_PARAM,STRING);
5926         para_name_list=get_list_of_tokens(para_name,NULL, &n_para_name);
5927         for ( a=0; a< n_para_name; a++)
5928             {
5929             sprintf ( parameter_list[number_of_parameters++],"%s", para_name_list[a]);
5930             }
5931         free_char(para_name_list,-1);
5932
5933
5934
5935
5936
5937         set_flag[0]=0;
5938         va_start (ap, usage);
5939
5940         if (strm3 (type, "S","R_F","W_F"))
5941                 string_val=va_arg(ap, char**);
5942         else if (strm2 (type, "D","FL"))
5943                 int_val=va_arg(ap, int*);
5944         else if (strm (type, "F"))
5945                 float_val=va_arg(ap, float*);
5946         else
5947             myexit (EXIT_FAILURE);
5948
5949
5950
5951         default_value1=va_arg(ap, char*);
5952         default_value2=va_arg(ap, char*);
5953         range_left    =va_arg(ap, char*);
5954         range_right   =va_arg(ap, char*);
5955         va_end(ap);
5956
5957
5958         para_name_list=get_list_of_tokens(para_name, NULL, &n_para_name);
5959         for ( a=0; a<n_para_name; a++)
5960             {
5961             if ( (pos=name_is_in_list(para_name_list[a], argv,argc, STRING))!=-1)break;
5962             }
5963         free_char (para_name_list,-1);
5964
5965         if (  (name_is_in_list("-help" , argv,argc  ,STRING)!=-1) && (argc==2 || (name_is_in_list( para_name , argv,argc  ,STRING)!=-1)))
5966           {
5967
5968                fprintf ( stderr, "PARAMETER   : %s\n",  para_name);
5969                fprintf ( stderr, "USAGE       : %s\n",       usage);
5970                fprintf ( stderr, "MAX_N_VALUES: %d\n",   max_n_val);
5971                fprintf ( stderr, "DEFAULT     : %s OR %s (when flag set)\n", default_value1, default_value2);
5972                fprintf ( stderr, "RANGE       : [%s]...[%s]\n", range_left,(strm(range_right,"any"))?"":range_right);
5973                fprintf ( stderr, "TYPE        : %s\n\n", type);
5974                return 0;
5975           }
5976         else if ( name_is_in_list ("-help" , argv,argc  ,STRING)!=-1)
5977           {
5978             return 0;
5979           }
5980         else if (para_name[0]!='-')
5981            {
5982            fprintf ( stderr, "\nWRONG PARAMETER DEFINITION %s Must Start with a dash", para_name);
5983            myexit (EXIT_FAILURE);
5984            }
5985         else if (pos==-1)
5986            {
5987            if ( optional==OPTIONAL)
5988               {
5989               set_flag[0]=0;
5990               para_val=get_list_of_tokens(default_value1, NULL, &n_para_val);
5991
5992               for (n_para=0; n_para<n_para_val && !strm (default_value1, "NULL"); n_para++)
5993                   {
5994                   if ( strm (para_val[n_para], ""))
5995                       {
5996                       set_flag[0]=0;
5997                       break;
5998                       }
5999                   else if ( strm (type, "FL"))
6000                       {
6001                       set_flag[0]=atoi(para_val[n_para]);
6002                       break;
6003                       }
6004                   else if (strm3 (type, "S", "R_F","W_F"))
6005                       {
6006
6007                       sprintf ( string_val[n_para], "%s",para_val[n_para]);
6008                       }
6009                   else if ( strm (type, "D"))
6010                       int_val[n_para]=atoi(para_val[n_para]);
6011                   else if ( strm (type, "F"))
6012                       float_val[n_para]=atof(para_val[n_para]);
6013                   }
6014               free_char (para_val, -1);
6015
6016               if (n_para==0 && strm3(type, "S","W_F","R_F") && strm (default_value1, "NULL"))
6017                   {
6018                   vfree (string_val[0]);
6019                   string_val[0]=NULL;
6020
6021                   }
6022               else if (n_para==0 && strm (type, "D") && strm (default_value1, "NULL"))int_val[0]=0;
6023               else if (n_para==0 && strm (type, "F") && strm (default_value1, "NULL"))float_val[0]=0;
6024
6025               }
6026            else
6027               {
6028               fprintf ( stderr, "\nParameter %s is not optional",para_name);
6029               myexit (EXIT_FAILURE);
6030               }
6031            }
6032         else if (pos!=-1)
6033           {
6034           set_flag[0]=1;
6035           for (a=pos+1; a< argc; a++)
6036               {
6037               if ( is_parameter(argv[a]))break;
6038               else
6039                   {
6040                   if ( n_para>=max_n_val)
6041                      {
6042                      n_para=max_n_val-1;
6043
6044                      }
6045                   if ( !(strm ( argv[a], "NULL")))
6046                     {
6047                       if ( strm3(type, "S", "R_F", "W_F"))
6048                           {
6049                           sprintf ( string_val[n_para],"%s", argv[a]);
6050                           }
6051                        else if (strm (type, "D"))
6052                           {
6053                           int_val[n_para]=atoi(argv[a]);
6054                           }
6055                        else if (strm ( type,"F"))
6056                           {
6057                           float_val[n_para]=atof(argv[a]);
6058                           }
6059                     }
6060                   n_para++;
6061                   }
6062               }
6063
6064           if ( n_para==0 && !strm2(default_value2,"","NULL") && !strm(type, "FL"))
6065               {
6066               para_val=get_list_of_tokens(default_value2, NULL, &n_para_val);
6067               for ( n_para=0; n_para<n_para_val; n_para++)
6068                   {
6069                   if ( strm3(type, "S", "R_F", "W_F"))sprintf ( string_val[n_para],"%s", para_val[n_para]);
6070                   else if (strm (type, "D"))int_val  [n_para]=atoi(para_val[n_para]);
6071                   else if (strm ( type,"F"))float_val[n_para]=atof(para_val[n_para]);
6072                   }
6073               free_char (para_val,-1);
6074               }
6075           else if (n_para==0 && strm (type, "FL"));
6076           else if (n_para==0 && strm3(type, "S","W_F","R_F") && strm (default_value2, "NULL")){vfree (string_val[0]);string_val[0]=NULL;}
6077           else if (n_para==0 && strm (type, "D") && strm (default_value2, "NULL"))int_val[0]=0;
6078           else if (n_para==0 && strm (type, "F") && strm (default_value2, "NULL"))float_val[0]=0;
6079           else if (n_para==0 && strm (default_value2, ""))
6080                  {
6081                  fprintf ( stderr, "\nParam %s needs a value [FATAL/%s]", para_name, argv[0]);
6082                  myexit(EXIT_FAILURE);
6083                  }
6084           else;
6085           }
6086
6087 /*Check That The Parameters are in the Good Range*/
6088
6089         pv_l=get_list_of_tokens( range_left , NULL, &n_pv_l);
6090         pv_r=get_list_of_tokens( range_right, NULL, &n_pv_r);
6091
6092         for ( a=0; a< n_para; a++)
6093             {
6094             if ( strm (type, "R_F") && !check_file_exists(string_val[a]) && !check_file_exists(string_val[a]+1))
6095                         {
6096                         fprintf ( stderr, "PARAM %s: File %s does not exist [FATAL/%s]\n",para_name,string_val[a], argv[0]);
6097                         myexit(EXIT_FAILURE);
6098                         }
6099             else if ( strm (pv_l[0], "any"));
6100             else if ( strm (type, "D"))
6101                  {
6102                  if ( n_pv_l==1)
6103                     {
6104                     min=(double)atoi(pv_l[0]);
6105                     max=(double)atoi(pv_r[0]);
6106                     if ( int_val[a]<min || int_val[a]>max)
6107                        {
6108                        fprintf ( stderr, "\n%s out of range [%d %d] [FATAL/%s]\n", para_name, (int)min, (int)max,argv[0]);
6109                        myexit (EXIT_FAILURE);
6110                        }
6111                     }
6112                  else
6113                     {
6114                     sprintf ( value, "%d", int_val[a]);
6115                     if ( name_is_in_list(value, pv_l, n_pv_l, STRING)==-1)
6116                         fprintf ( stderr, "\n%s out of range [%s: ", para_name, value);
6117                     print_array_char (stderr, pv_l, n_pv_l, " ");
6118                     fprintf ( stderr, "\n");
6119                     myexit(EXIT_FAILURE);
6120                     }
6121                  }
6122             else if ( strm (type, "F"))
6123                  {
6124                   if ( n_pv_l==1)
6125                     {
6126                     min=(double)atof(range_left);
6127                     max=(double)atof(range_right);
6128                     if ( float_val[a]<min || float_val[a]>max)
6129                        {
6130                        fprintf ( stderr, "\n%s out of range [%f %f] [FATAL/%s]\n", para_name, (float)min, (float)max,argv[0]);
6131                        myexit (EXIT_FAILURE);
6132                        }
6133                      }
6134                   else
6135                      {
6136                      sprintf ( value, "%f", float_val[a]);
6137                      if ( name_is_in_list(value, pv_l, n_pv_l, STRING)==-1)
6138                         fprintf ( stderr, "\n%s out of range [%s: ", para_name, value);
6139                      print_array_char (stderr, pv_l, n_pv_l, " ");
6140                      fprintf ( stderr, "\n");
6141                      myexit(EXIT_FAILURE);
6142                      }
6143                  }
6144             }
6145
6146
6147          if ( fp[0]!=NULL)
6148               {
6149               fprintf (fp[0], "%-15s\t%s\t[%d] ", para_name, type, set_flag[0]);
6150               for (a=0; a<n_para; a++)
6151                   {
6152                   if ( strm3 ( type, "S", "R_F", "W_F"))fprintf ( fp[0], "\t%s", string_val[a]);
6153                   else if ( strm  ( type, "D"))fprintf ( fp[0], "\t%d ", int_val[a]);
6154                   else if ( strm  ( type, "F"))fprintf ( fp[0], "\t%f ", float_val[a]);
6155                  }
6156               if ( strm (type, "FL"))fprintf ( fp[0], "\t%d", int_val[0]);
6157               fprintf ( fp[0], "\n");
6158               }
6159
6160         free_char ( pv_l, -1);
6161         free_char ( pv_r, -1);
6162         return n_para;
6163         }
6164
6165
6166
6167 char ** get_parameter ( char *para_name, int *np, char *fname)
6168 {
6169     /*
6170     In:
6171     para_name: the name of the parameter to look for
6172     fname: the name of the file containing the parameters
6173     np[0]: set to 0
6174
6175     Out:
6176     char ** containing the np[0] values taken by para_name in fname.
6177
6178     Special:
6179     if fname=NULL, para_name is searched using the last value taken by fp.
6180
6181     Note: by default, the function keeps a file handle open until the first unsuccessful call.
6182     */
6183
6184     static FILE *fp;
6185     static char *line;
6186     char ** return_value;
6187
6188     if ( strm (para_name, "CLOSE THE FILE"))
6189       {
6190         vfclose ( fp);
6191         return NULL;
6192       }
6193
6194     if ( line==NULL)line=vcalloc ( VERY_LONG_STRING+1, sizeof (char));
6195     if ( fname!=NULL && fp!=NULL)vfclose (fp);
6196
6197     np[0]=0;
6198
6199     if ((fp=find_token_in_file ( fname,(fname==NULL)?fp:NULL, para_name))==NULL)
6200         {
6201              return NULL;
6202         }
6203     else
6204         {
6205         fgets ( line, VERY_LONG_STRING,fp);
6206         return_value=get_list_of_tokens ( line, NULL, np);
6207         return return_value;
6208         }
6209 }
6210
6211 FILE * set_fp_id ( FILE *fp, char *id)
6212         {
6213 /*Sets fp just after id, id needs to be at the begining of the line*/
6214         char string[10000];
6215         int cont=1;
6216         int c;
6217
6218         while ( cont==1)
6219                 {
6220                 c=fgetc(fp);
6221                 if ( c!=EOF)
6222                         {
6223
6224                         ungetc(c, fp);
6225                         fscanf ( fp, "%s", string);
6226
6227                         if ( strcmp ( string, id)==0)
6228                                 return fp;
6229                         else while ( c!='\n' && c!=EOF)
6230                                 c=fgetc(fp);
6231                         }
6232                 else if ( c==EOF)
6233                         {
6234                         fclose ( fp);
6235                         return NULL;
6236                         }
6237                 }
6238         return fp;
6239         }
6240 FILE * set_fp_after_char ( FILE *fp, char x)
6241         {
6242 /*sets fp just after the first occurence of x*/
6243
6244
6245
6246         int cont=1;
6247         int c;
6248
6249         while ( cont==1)
6250                 {
6251                 c=fgetc(fp);
6252                 if ( c!=EOF)
6253                         {
6254                         if ( c==x)
6255                                 return fp;
6256
6257                         }
6258                 else if ( c==EOF)
6259                         {
6260                         fclose ( fp);
6261                         return NULL;
6262                         }
6263                 }
6264         return fp;
6265         }
6266
6267 FILE * find_token_in_file_nlines ( char *fname, FILE * fp, char *token, int n_line)
6268
6269         {
6270           /*This function finds the string TOKEN (as a single word) in the n_line first lines of fname.
6271             It returns NULL if not found or the position of the fp
6272           */
6273
6274           char *tmp_name=NULL;
6275           FILE *fp1;
6276           FILE *fp2;
6277           char buffer[10001];
6278           int a;
6279           if ( !fp && !check_file_exists(fname))return NULL;
6280           if ( fp==NULL)
6281             {
6282               tmp_name=vtmpnam ( NULL);
6283
6284               fp1=vfopen (fname, "r");
6285               fp2=vfopen (tmp_name, "w");
6286
6287               for ( a=0; a< n_line && fgets(buffer, 10000, fp1)!=NULL; a++)
6288                 {
6289                   fprintf ( fp2, "%s", buffer);
6290                 }
6291               vfclose (fp1);
6292               vfclose (fp2);
6293             }
6294           return find_token_in_file ( tmp_name,fp,token);
6295         }
6296
6297 int token_is_in_file ( char *fname, char *token)
6298 {
6299   /*TH:an argument against torture: innocents get tortured longer, likewise for token testing*/
6300   FILE *fp;
6301   static char *buf;
6302   char *b, *p;
6303   int begining;
6304
6305
6306
6307   if (token[0]=='\n')
6308     {
6309
6310       begining=1;
6311       token++;
6312     }
6313   else
6314     {
6315       begining=0;
6316     }
6317
6318   if ( !fname || !check_file_exists(fname))return 0;
6319   else
6320     {
6321       fp=vfopen (fname, "r");
6322       while ((b=vfgets (buf,fp))!=NULL)
6323         {
6324           buf=b;
6325           p=strstr (buf, token);
6326           if (!p);
6327           else if ( begining==1 && p==buf){vfclose (fp); return 1;}
6328           else if ( begining==0 && p){vfclose (fp); return 1;}
6329         }
6330       vfclose (fp);
6331       return 0;
6332     }
6333   return 0;
6334 }
6335
6336 char *vfgets ( char *buf, FILE *fp)
6337 {
6338   int c;
6339   int buf_len, l;
6340   char *bufin;
6341   static int debug;
6342
6343
6344   if ( !debug)
6345     debug=(getenv ("DEBUG_TCOFFEE"))?1:-1;
6346
6347
6348   bufin=buf;
6349
6350   if (!buf)
6351     {
6352       buf=vcalloc ( 1000, sizeof (char));
6353       buf_len=1000;
6354     }
6355   else
6356     {
6357       buf_len=read_array_size_new (buf);
6358     }
6359
6360   l=0;
6361
6362   if ( (c=fgetc (fp))==EOF)return NULL;
6363   else ungetc (c, fp);
6364
6365
6366   while ( (c=fgetc (fp))!='\n' && c!=EOF)
6367     {
6368       if (l>=buf_len)
6369         {buf_len+=100;buf=vrealloc (buf, buf_len*sizeof (char));}
6370       buf[l++]=c;
6371     }
6372   /*Add the cariage return*/
6373   if ( c=='\n')
6374     {
6375       if (l>=buf_len){buf_len+=100,buf=vrealloc (buf, buf_len*sizeof (char));}
6376       buf[l++]='\n';
6377     }
6378   /*add the terminator*/
6379   if (l>=buf_len){buf_len+=100,buf=vrealloc (buf, buf_len*sizeof (char));}
6380   buf[l]='\0';
6381
6382   if ( bufin!=buf && bufin!=NULL && debug==1)
6383     fprintf ( stderr, "\nPointer change in vfgets...");
6384
6385   return buf;
6386 }
6387
6388
6389 FILE * find_token_in_file ( char *fname, FILE * fp, char *token)
6390         {
6391         int c;
6392         static char *name;
6393         int token_len;
6394
6395         int only_start;
6396
6397         /*Note: Token: any string
6398                 If Token[0]=='\n' Then Token only from the beginning of the line
6399         */
6400
6401         if (!fp && !check_file_exists(fname))return NULL;
6402
6403         if ( token[0]=='\n'){token++;only_start=1;}
6404         else only_start=0;
6405
6406         token_len=strlen (token);
6407
6408
6409
6410
6411
6412         if (!fp)
6413           {
6414             if (name)vfree (name);
6415             name = vcalloc (((fname)?measure_longest_line_in_file (fname):10000)+1, sizeof (char));
6416             fp=vfopen ( fname, "r");
6417           }
6418
6419         while ( (fscanf ( fp, "%s", name))!=EOF)
6420                 {
6421
6422                 if ( name[0]=='*')while ( ((c=fgetc (fp))!='\n')&& (c!=EOF));
6423                 else if (strncmp ( name, token,token_len)==0){return fp;}
6424                 else if (only_start) while ( ((c=fgetc (fp))!='\n')&& (c!=EOF));
6425                 }
6426
6427         vfclose ( fp);
6428         return NULL;
6429         }
6430 int **get_file_block_pattern (char *fname, int *n_blocks, int max_n_line)
6431         {
6432         int c;
6433         FILE *fp;
6434         char *line;
6435         int lline;
6436         int **l;
6437         int in_block;
6438
6439         int max_block_size;
6440         int block_size;
6441         int x;
6442         int n_line;
6443
6444         lline=measure_longest_line_in_file (fname)+1;
6445         line=vcalloc ( sizeof (char),lline+1);
6446
6447         fp=vfopen (fname, "r");
6448         max_block_size=block_size=0;
6449         in_block=1;
6450         n_blocks[0]=0;
6451         n_line=0;
6452         while ((c=fgetc(fp))!=EOF && (n_line<max_n_line || !max_n_line))
6453             {
6454                   ungetc (c, fp);
6455                   fgets ( line, lline,fp);
6456                   n_line++;
6457
6458                   if ( is_alnum_line (line) && !in_block){n_blocks[0]++;in_block=1;}
6459                   if ( is_alnum_line (line))
6460                      {
6461                      block_size++;
6462                      }
6463                   else
6464                       {
6465                       in_block=0;
6466                       max_block_size=MAX( max_block_size, block_size);
6467                       block_size=0;
6468                       }
6469             }
6470
6471
6472         max_block_size=MAX( max_block_size, block_size);
6473         vfclose ( fp);
6474
6475         l=declare_int (n_blocks[0]+1,max_block_size+1);
6476
6477
6478         fp=vfopen (fname, "r");
6479         in_block=1;
6480         n_blocks[0]=0;
6481         n_line=0;
6482         while ((c=fgetc(fp))!=EOF && (n_line<max_n_line || !(max_n_line)))
6483               {
6484                   ungetc (c, fp);
6485                   fgets ( line, lline,fp);
6486                   n_line++;
6487
6488                   if ( is_alnum_line (line) && !in_block){n_blocks[0]++;in_block=1;}
6489                   if ( is_alnum_line (line))
6490                       {
6491                       l[n_blocks[0]][0]++;
6492                       free_char (get_list_of_tokens (line, " \t\n*:,", &x), -1);
6493
6494                       if ( l[n_blocks[0]][0]> max_block_size)fprintf ( stderr, "\nERROR %d", l[n_blocks[0]][0]);
6495
6496                       l[n_blocks[0]] [l[n_blocks[0]][0]]=x;
6497                       }
6498                   else
6499                       {
6500                           in_block=0;
6501                       }
6502               }
6503         n_blocks[0]++;
6504         vfree(line);
6505         vfclose (fp);
6506         return l;
6507         }
6508
6509 char * strip_file_from_comments (char *com, char *in_file)
6510 {
6511   /*Removes in file in_file every portion of line to the right of one of the symbols included in com
6512     Writes the striped file into a vtmpnam file
6513    */
6514   FILE *fp1;
6515   FILE *fp2;
6516   char *out_file;
6517   int c;
6518
6519   out_file=vtmpnam(NULL);
6520
6521
6522   fp1=vfopen (in_file , "r");
6523   fp2=vfopen (out_file, "w");
6524   while ( (c=fgetc(fp1))!=EOF)
6525           {
6526             if (strchr(com, c))
6527               {
6528                 while ( (c=fgetc(fp1))!='\n' && c!=EOF);
6529               }
6530             else
6531               {
6532                 fprintf (fp2, "%c", c);
6533                 while ( (c=fgetc(fp1))!='\n' && c!=EOF)fprintf (fp2, "%c", c);
6534                 if ( c!=EOF)fprintf (fp2, "%c", c);
6535               }
6536           }
6537   vfclose (fp1);
6538   vfclose (fp2);
6539
6540   return out_file;
6541 }
6542 FILE * skip_commentary_line_in_file ( char com, FILE *fp)
6543 {
6544   int c=0;
6545
6546   if ( fp==NULL)return NULL;
6547   while ((c=fgetc(fp))==com)
6548     {
6549       while ((c=fgetc(fp))!='\n' && c!=EOF);
6550     }
6551   if ( c!=EOF && c!='\n')ungetc(c, fp);
6552   return fp;
6553 }
6554
6555 int check_for_update ( char *web_address)
6556 {
6557   char command[1000];
6558   char *file;
6559   float new_version, old_version;
6560   FILE *fp;
6561
6562   check_internet_connection (IS_NOT_FATAL);
6563   file=vtmpnam(NULL);
6564
6565   sprintf ( command, "%s/%s.version",DISTRIBUTION_ADDRESS, PROGRAM);
6566   url2file ( command, file);
6567
6568   fp=vfopen ( file, "r");
6569   fscanf ( fp, "Version_%f", &new_version);
6570   vfclose ( fp);
6571   sscanf ( VERSION, "Version_%f", &old_version);
6572
6573   if ( old_version<new_version)
6574     {
6575       fprintf ( stdout, "\nUpdate Status: outdated");
6576       fprintf ( stdout, "\nYour version of %s is not up to date", PROGRAM);
6577       fprintf ( stdout, "\nDownload the latest version %.2f with the following command:", new_version);
6578       fprintf ( stdout, "\n\twget %s/%s_distribution_Version_%.2f.tar.gz\n", DISTRIBUTION_ADDRESS, PROGRAM, new_version);
6579       return EXIT_FAILURE;
6580     }
6581   else if ( old_version> new_version)
6582     {
6583       fprintf ( stdout, "\nUpdate Status: beta-release");
6584       fprintf ( stdout, "\nYour are using a beta-release of %s(%s)\n", PROGRAM, VERSION);
6585     }
6586   else
6587     {
6588       fprintf (stdout, "\nUpdate Status: uptodate");
6589       fprintf (stdout, "\nProgram %s(%s) is up to date\n", PROGRAM, VERSION);
6590     }
6591   return EXIT_SUCCESS;
6592 }
6593
6594
6595
6596
6597
6598 int check_environement_variable_is_set ( char *variable, char *description, int fatal)
6599 {
6600   if ( getenv (variable)==NULL)
6601     {
6602       fprintf ( stderr, "\nERROR: You must set %s\n%s %s", variable, description, description);
6603       if ( fatal==IS_FATAL)
6604         {
6605         fprintf ( stderr, "\n[%s:FATAL]\n", PROGRAM);
6606         exit (EXIT_FAILURE);
6607         }
6608       else
6609         fprintf ( stderr, "\n[%s:WARNING]\n", PROGRAM);
6610     }
6611   return 1;
6612 }
6613
6614 int url2file (char *address, char *out)
6615 {
6616   char command[1000];
6617
6618
6619   if      (check_program_is_installed ("wget",NULL, NULL,WGET_ADDRESS, IS_NOT_FATAL))sprintf (command, "wget %s -O%s >/dev/null 2>/dev/null", address, out);
6620   else if (check_program_is_installed ("curl",NULL, NULL,CURL_ADDRESS, IS_NOT_FATAL))sprintf (command, "curl %s -o%s >/dev/null 2>/dev/null", address, out);
6621   else
6622     {
6623       printf_exit (EXIT_FAILURE, stderr, "ERROR: Impossible to fectch external file: Neither wget nor curl is installed on your system [FATAL:%s]\n", PROGRAM);
6624       return EXIT_FAILURE;
6625     }
6626   return safe_system (command);
6627 }
6628
6629 int wget (char *address, char *out)
6630 {
6631   return printf_system ( "curl %s -O%s >/dev/null 2>/dev/null", address, out);
6632   }
6633
6634 int curl (char *address, char *out)
6635 {
6636   return printf_system ( "curl %s -o%s >/dev/null 2>/dev/null", address, out);
6637 }
6638
6639
6640 int simple_check_internet_connection (char *ref_site)
6641 {
6642   char *test,command[1000];
6643   int n, internet=0;
6644
6645   test=vtmpnam (NULL);
6646   if (url2file((ref_site)?ref_site:TEST_WWWSITE_4_TCOFFEE,test)!=EXIT_SUCCESS)internet=0;
6647   else if ((n=count_n_char_in_file(test))<10)internet=0;
6648   else internet =1;
6649
6650   return internet;
6651 }
6652 int check_internet_connection  (int mode)
6653 {
6654   int internet;
6655   internet=simple_check_internet_connection (NULL);
6656   if (internet)return 1;
6657   else if ( mode == NON_INTERACTIVE)return internet;
6658   else if ( mode == IS_FATAL)
6659     {
6660       add_warning ( stderr,"\nERROR: You do not seem to have an active Internet Connection. Check your Proxy Setting [proxy:%s][%s:SERIOUS]\n", getenv ("http_proxy"), PROGRAM);
6661       return EXIT_FAILURE;
6662     }
6663   else if ( mode==IS_NOT_FATAL)
6664     {
6665
6666       fprintf ( stderr, "\n\n\n");
6667       fprintf ( stderr, "*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
6668       fprintf ( stderr, "*\n");
6669       fprintf ( stderr, "*   Your Internet Connection Does not seem to be available. You may need to reconfigure your Proxy\n");
6670       fprintf ( stderr, "*   ------- Current Proxy Value: [%s]\n*\n", (getenv ("http_proxy")?getenv ("http_proxy"):"NO Proxy Set"));
6671       Proxy(INPUT, RESET);
6672
6673       return check_internet_connection (IS_FATAL);
6674     }
6675   else
6676     {
6677       return 0;
6678     }
6679 }
6680 char *pg2path (char *pg)
6681 {
6682   char *path;
6683   char *p;
6684   char *tmp;
6685   FILE *fp;
6686
6687   if ( !pg) return NULL;
6688   tmp=vtmpnam(NULL);
6689
6690   printf_system_direct("which %s>%s 2>/dev/null", pg, tmp);
6691   path=file2string (tmp);
6692   chomp (path);
6693   if (!check_file_exists (path) && !strstr (pg, ".exe"))
6694     {
6695       char pg2[1000];
6696       sprintf ( pg2, "%s.exe", pg);
6697       path=pg2path (pg2);
6698     }
6699
6700   return path;
6701 }
6702
6703
6704 int check_program_is_installed ( char *program_name, char *path_variable, char *path_variable_name, char *where2getit, int fatal)
6705   {
6706    char command[LONG_STRING];
6707    static char *path;
6708
6709    if ( strm (where2getit, "built_in"))return 1;
6710
6711    if (path)vfree (path);
6712
6713    if ( check_file_exists (path_variable))
6714      {
6715        return 1;
6716      }
6717    else
6718      {
6719
6720        path=pg2path (program_name);
6721        if (path && path[0])return 1;
6722        else
6723          {
6724            int install=EXIT_FAILURE;
6725
6726            if (fatal==INSTALL || fatal==INSTALL_OR_DIE)
6727              {
6728                HERE ("************** %s is missing from your system. T-Coffee will make an attempt to install it.\n", program_name);
6729                install=printf_system ("install.pl %s -plugins=%s -clean", program_name, get_plugins_4_tcoffee(NULL));
6730              }
6731            if ( install==EXIT_SUCCESS)return 1;
6732            else if ( fatal==INSTALL)return 0;
6733            else if ( fatal==NO_REPORT)return 0;
6734
6735            if (fatal==IS_FATAL || fatal==INSTALL_OR_DIE)check_configuration4program();
6736
6737            fprintf ( stderr, "\n#*****************************************************************");
6738            if (fatal) fprintf ( stderr, "\n#ERROR [FATAL:%s]", PROGRAM);
6739            else fprintf ( stderr, "\n#WARNING [%s]", PROGRAM);
6740            fprintf ( stderr, "\n# The Program %s Needed by %s Could not be found", program_name, PROGRAM);
6741            fprintf ( stderr, "\n# If %s is installed on your system:", program_name);
6742            fprintf ( stderr, "\n#\t     -Make sure %s is in your $path:",program_name);
6743
6744            fprintf ( stderr, "\n# If %s is NOT installed obtain a copy from:", program_name);
6745            fprintf ( stderr, "\n#\t%s\n#\n#",where2getit);
6746            fprintf ( stderr, "\n# and install it manualy");
6747            fprintf ( stderr, "\n******************************************************************\n");
6748          }
6749      }
6750    if ( fatal==IS_FATAL || fatal==INSTALL_OR_DIE) myexit (EXIT_FAILURE);
6751    return 0;
6752   }
6753
6754 FILE * display_output_filename ( FILE *io, char *type, char *format, char *name, int check_output)
6755 {
6756   static char ***buf;
6757   static int nbuf;
6758
6759   if ( strm ( name, "stdout") || strm (name, "stderr"))return io;
6760
6761   if ( check_output==STORE)
6762     {
6763       int a;
6764       if ( buf==NULL)buf=vcalloc ( 1000, sizeof (char**));
6765
6766       for (a=0; a<nbuf; a++)
6767         if ( strm (name, buf[a][2]))return io;
6768
6769       buf[nbuf]=declare_char (3, 1000);
6770       sprintf ( buf[nbuf][0], "%s", type);
6771       sprintf ( buf[nbuf][1], "%s", format);
6772       sprintf ( buf[nbuf][2], "%s", name);
6773       nbuf++;
6774       return io;
6775     }
6776   else if ( check_output==FLUSH)
6777     {
6778       int a;
6779       for ( a=0; a< nbuf; a++)
6780         {
6781           io=display_output_filename ( io, buf[a][0], buf[a][1], buf[a][2], CHECK);
6782           free_char (buf[a], -1);
6783         }
6784       nbuf=0;
6785     }
6786   else if ( check_output==CHECK)
6787     {
6788       if (check_file_exists(name)==NULL)
6789         {
6790           fprintf ( io, "\n\t#### File Type= %10s Format= %10s Name= %s | NOT PRODUCED [WARNING:%s:%s]\n",type, format, name, PROGRAM, VERSION );
6791           return io;
6792         }
6793       add_file2error_report (type, format, name);
6794       fprintf ( io, "\n\t#### File Type= %10s Format= %10s Name= %s",type, format, name );
6795     }
6796   return io;
6797 }
6798
6799 FILE * display_input_filename ( FILE *io, char *type, char *format, char *name, int check_output)
6800 {
6801   if ( check_output==CHECK && check_file_exists(name)==NULL)
6802     {
6803       fprintf ( io, "\n\tIIII INPUT File Type= %10s Format= %10s Name= %s | NOT PRODUCED [WARNING:%s:%s]\n",type, format, name, PROGRAM, VERSION );
6804       return io;
6805     }
6806   add_file2error_report (type, format, name);
6807
6808   fprintf ( io, "\n\t#### File Type= %10s Format= %10s Name= %s",type, format, name );
6809   return io;
6810 }
6811
6812 int file_is_empty (char *fname)
6813 {
6814   struct stat s;
6815   if (!fname) return 1;
6816
6817   stat (fname, &s);
6818   if (s.st_size)return 0;
6819   else return 1;
6820   }
6821 int file_exists (char *path, char *fname)
6822 {
6823   struct stat s;
6824   char file[1000];
6825
6826   if (!fname)return 0;
6827   else if (!path)sprintf (file, "%s", fname);
6828   else sprintf ( file, "%s/%s", path, fname);
6829
6830
6831
6832   if (stat(file,& s)!=-1)
6833     return S_ISREG(s.st_mode);
6834   else return 0;
6835 }
6836
6837 int isdir  (char *file)
6838 {
6839   struct stat s;
6840   if (stat (file,&s)!=-1)
6841     return S_ISDIR(s.st_mode);
6842   else return 0;
6843 }
6844 int rrmdir (char *s)
6845 {
6846   if (isdir(s))return printf_system_direct ("rm -r %s", s);
6847   return EXIT_FAILURE;
6848 }
6849
6850 int isexec (char *file)
6851 {
6852   char *state;
6853
6854   state=ls_l(NULL,file);
6855   if (state[0]==0) return 0;
6856
6857   if ( state[0]=='d') return 0;
6858   if ( state[3]=='x') return 1;
6859   if ( state[6]=='x') return 1;
6860   if ( state[9]=='x') return 1;
6861   return 0;
6862 }
6863
6864 char *ls_l ( char *path,char *file)
6865 {
6866   char *tmpfile;
6867   static char *state;
6868   char command[10000];
6869   FILE *fp;
6870   int a;
6871
6872   tmpfile=vtmpnam (NULL);
6873   if (!state)
6874     {
6875
6876       state=vcalloc (100, sizeof (char));
6877     }
6878   for (a=0;a<100; a++)state[a]=0;
6879   if (!file)return state;
6880
6881
6882   sprintf ( command, "ls -l %s%s%s >%s 2>/dev/null", (path!=NULL)?path:"", (path!=NULL)?"/":"",file, tmpfile);
6883   safe_system (command);
6884   fp=vfopen (tmpfile, "r");
6885   if (!fscanf ( fp, "%s", state))
6886     {
6887       vfclose(fp); return 0;
6888     }
6889   vfclose (fp);
6890   return state;
6891 }
6892
6893 int my_mkdir ( char *dir_in)
6894 {
6895
6896   int dir_sep='/';
6897
6898   int a, buf;
6899   char *dir;
6900
6901   dir=vcalloc ( strlen (dir_in)+strlen (get_home_4_tcoffee())+100, sizeof (char));
6902   sprintf ( dir, "%s", dir_in);
6903   tild_substitute ( dir, "~",get_home_4_tcoffee());
6904
6905
6906
6907   a=0;
6908
6909   while (dir[a]!='\0')
6910     {
6911
6912       if ( dir[a]==dir_sep || dir[a+1]=='\0')
6913         {
6914           buf= dir[a+1];
6915           dir[a+1]='\0';
6916
6917           if (access(dir, F_OK)==-1)
6918             {
6919
6920               printf_system_direct("mkdir %s", dir);
6921
6922               if ( access (dir, F_OK)==-1)
6923                 {
6924                   fprintf ( stderr, "\nERROR: Could Not Create Directory %s [FATAL:%s]", dir, PROGRAM);
6925                   exit (EXIT_FAILURE);
6926                 }
6927             }
6928           dir[a+1]=buf;
6929         }
6930       a++;
6931     }
6932
6933   vfree (dir);
6934   return 1;
6935 }
6936
6937 int filename_is_special (char *fname)
6938 {
6939   if ( strm5 (fname, "default", "stdin", "stdout","stderr", "/dev/null"))return 1;
6940   if ( strm3 (fname, "STDIN", "STDOUT", "STDERR"))return 1;
6941   return 0;
6942 }
6943
6944 char* check_file_exists ( char *fname_in)
6945         {
6946         FILE *fp;
6947         static char *fname1;
6948         static char *fname2;
6949
6950
6951         if (!fname_in)return NULL;
6952         if (!fname_in[0])return NULL;
6953         if (fname_in[0]=='-')return NULL;
6954
6955         if (!fname1){fname1=vcalloc (1000, sizeof (char));}
6956         if (!fname2){fname2=vcalloc (1000, sizeof (char));}
6957
6958         sprintf ( fname1, "%s", fname_in);tild_substitute (fname1, "~", get_home_4_tcoffee());
6959         sprintf ( fname2, "%s%s", get_cache_dir(),fname1);
6960
6961         if ( filename_is_special (fname1))return fname1;
6962         if ( strm5 (fname1, "no", "NO", "No", "NO_FILE","no_file"))return NULL/*fname1*/;
6963         if (!file_exists( NULL,fname1))
6964           {
6965             if (!file_exists (NULL,fname2))return NULL;
6966             else return fname2;
6967           }
6968         else return fname1;
6969         return NULL;
6970         }
6971
6972
6973 void create_file ( char *name)
6974         {
6975         FILE *fp;
6976
6977         fp=fopen (name, "w");
6978         fclose (fp);
6979         }
6980 void delete_file ( char *fname)
6981         {
6982         char command[1000];
6983         FILE * fp;
6984
6985         fp=fopen ( fname, "w");
6986         fprintf ( fp, "x");
6987         fclose ( fp);
6988
6989         sprintf ( command, "rm %s", fname);
6990         my_system ( command);
6991
6992         }
6993
6994 int util_rename ( char *from, char *to)
6995         {
6996         FILE *fp_from;
6997         FILE *fp_to;
6998         int c;
6999
7000
7001         if ( !check_file_exists (from))return 0;
7002         else if ( check_file_exists (to) && !vremove (to) && !rename ( from, to)==0 );
7003         else
7004                 {
7005
7006                 fp_from=vfopen ( from, "r");
7007                 fp_to=vfopen ( to, "w");
7008
7009                 while ( (c=fgetc (fp_from))!=EOF)fprintf ( fp_to, "%c", c);
7010
7011                 fclose (fp_from);
7012                 fclose ( fp_to);
7013
7014                 vremove ( from);
7015                 return 1;
7016                 }
7017         return 0;
7018         }
7019
7020
7021 int util_copy (  char *from, char *to)
7022         {
7023         FILE *fp_from;
7024         FILE *fp_to;
7025         int c;
7026
7027
7028         if (!check_file_exists (from))return 0;
7029         else
7030                 {
7031
7032                 fp_from=vfopen ( from, "r");
7033                 fp_to=vfopen ( to, "w");
7034
7035                 while ( (c=fgetc (fp_from))!=EOF)fprintf ( fp_to, "%c", c);
7036
7037                 fclose (fp_from);
7038                 fclose ( fp_to);
7039                 return 1;
7040                 }
7041         return 0;
7042         }
7043 FILE * output_completion4halfmat ( FILE *fp,int n, int tot, int n_reports, char *s)
7044
7045 {
7046   int max, left, achieved;
7047   int up;
7048
7049   if (n>=0)up=1;
7050   else up=-1;
7051
7052
7053   max=((tot*tot)-tot)/2;
7054   left=((tot-n)*(tot-n)-(tot-n))/2;
7055
7056   achieved=max-left;
7057   if (up==1);
7058   else
7059     {
7060       int b;
7061       b=achieved;
7062       achieved=left;
7063       left=b;
7064     }
7065   return output_completion (fp,achieved, max, n_reports, s);
7066 }
7067
7068
7069 FILE * output_completion ( FILE *fp,int n, int tot, int n_reports, char *string)
7070         {
7071
7072           static int ref_val;
7073           static int flag;
7074           static int ref_time;
7075           int t, elapsed;
7076           n++;
7077
7078           if ( n==1)
7079             {
7080               ref_val=flag=0;
7081               ref_time=get_time()/1000;
7082             }
7083           t=get_time()/1000;
7084           elapsed=t-ref_time;
7085
7086           if ( !ref_val && !flag)
7087             {
7088               fprintf (fp, "\n\t\t[%s][TOT=%5d][%3d %%][ELAPSED TIME: %4d sec.]",(string)?string:"",tot,(tot==1)?100:0, elapsed);
7089               flag=1;
7090             }
7091           else if ( n==tot)fprintf (fp, "\r\t\t[%s][TOT=%5d][%3d %%][ELAPSED TIME: %4d sec.]",(string)?string:"", tot,100, elapsed);
7092           else if ( ((n*100)/tot)>ref_val)
7093             {
7094               ref_val=((n*100)/tot);
7095               t=(ref_val==0)?0:elapsed/ref_val;
7096               t=t*(100-ref_val);
7097               t=0;
7098               fprintf (fp, "\r\t\t[%s][TOT=%5d][%3d %%][ELAPSED TIME: %4d sec.]", (string)?string:"",tot,ref_val, elapsed);
7099               flag=0;
7100             }
7101           return fp;
7102         }
7103 void * null_function (int a,...)
7104 {
7105   fprintf ( stderr, "\n[ERROR] Attempt to use the Null Function [FATAL:%s]", PROGRAM);
7106   crash ("");
7107   return NULL;
7108 }
7109
7110 int  btoi ( int nc,...)
7111 {
7112   va_list ap;
7113   int a, b;
7114   va_start (ap, nc);
7115   for ( a=0, b=0; a< nc; a++)
7116     {
7117       b+=pow(2,a)*va_arg (ap,int);
7118     }
7119   va_end(ap);
7120   return b;
7121 }
7122
7123 /*********************************************************************/
7124 /*                                                                   */
7125 /*                         Geometric FUNCTIONS                    */
7126 /*                                                                   */
7127 /*                                                                   */
7128 /*********************************************************************/
7129
7130 float get_geometric_distance ( float ** matrix, int ncoor, int d1, int d2, char *mode)
7131 {
7132   float d;
7133   float t=0;
7134   int a;
7135
7136   if ( strm (mode, "euclidian"))
7137     {
7138       for ( a=0; a< ncoor; a++)
7139         {
7140           d=(matrix[d1][a]-matrix[d2][a]);
7141           t+=d*d;
7142         }
7143       return (float)sqrt((double)t);
7144     }
7145   return 0;
7146 }
7147
7148
7149
7150 /*********************************************************************/
7151 /*                                                                   */
7152 /*                         MATHEMATICAL FUNCTIONS                    */
7153 /*                                                                   */
7154 /*                                                                   */
7155 /*********************************************************************/
7156 static double EXP_UNDERFLOW_THRESHOLD = -4.60f;
7157 static double LOG_UNDERFLOW_THRESHOLD = 7.50f;
7158 static double LOG_ZERO = -FLT_MAX;
7159 static double LOG_ONE = 0.0f;
7160 double log_addN (int N, double*L)
7161
7162 {
7163   double v;
7164   int a;
7165   if (N==0)return 0;
7166   if ( N==1)return L[0];
7167
7168   v=L[0];
7169   for ( a=1; a<N; a++)
7170     {
7171       v=log_add2(v, L[a]);
7172     }
7173   return v;
7174   }
7175 double log_add6 (double x, double y, double z, double w, double v, double e)
7176 {
7177   return log_add2(log_add3(x, y, z),log_add3(z,w,e));
7178 }
7179 double log_add5 (double x, double y, double z, double w, double v)
7180 {
7181   return log_add3(log_add2(x, y),log_add2(z,w),v );
7182 }
7183 double log_add4 (double x, double y, double z, double w)
7184 {
7185   return log_add2(log_add2(x, y),log_add2(z,w));
7186 }
7187 double log_add3 (double x, double y, double z)
7188 {
7189   return log_add2(log_add2(x, y),z);
7190 }
7191 double log_add2 ( double x, double y)
7192 {
7193   if (x < y)
7194     x = (x == LOG_ZERO || ((y - x) >= LOG_UNDERFLOW_THRESHOLD)) ? y : log (exp (x-y) + 1) + x;
7195   else
7196     x = (y == LOG_ZERO || ((x - y) >= LOG_UNDERFLOW_THRESHOLD)) ? x : log (exp (x-y) + 1) + y;
7197   return x;
7198 }
7199
7200
7201
7202
7203 float M_chooses_Nlog ( int m, int N)
7204 {
7205   /*Choose M elemets in N*/
7206   float  z1, z2,z=0;
7207   if ( m==N) return 0;
7208   else if ( m>N)
7209     {
7210       fprintf ( stderr, "\nERROR: M chosses N out of bounds ( M>N) [FATAL:%s]", PROGRAM);
7211       myexit (EXIT_FAILURE);
7212     }
7213   else
7214     {
7215       z1=factorial_log (m+1, N);
7216       z2=factorial_log (1, N-m);
7217       z=z1-z2;
7218       return z;
7219     }
7220
7221   return -1;
7222 }
7223
7224 float factorial_log ( int start, int end)
7225 {
7226   if ( end==0)return 0;
7227   else if ( end==start) return (float)my_int_log((double)start);
7228   else if ( start>end)
7229      {
7230        fprintf ( stderr, "\nERROR: factorial log out of bounds (%d %d) [FATAL:%s]",start, end, PROGRAM);
7231        myexit (EXIT_FAILURE);
7232     }
7233   else
7234     {
7235       int a=0;
7236       float x=0;
7237       for ( x=0,a=start; a<=end; a++)
7238         {
7239           x+=(float)my_int_log(a);
7240         }
7241       return x;
7242     }
7243   return 0;
7244 }
7245
7246 float my_int_log(int a)
7247 {
7248
7249   if ( a>=100000)return log(a);
7250   else
7251     {
7252       static float *lu;
7253       if (!lu) lu=vcalloc ( 100000, sizeof (float));
7254       if ( !lu[a]){lu[a]=log(a);}
7255       return lu[a];
7256     }
7257   return 0;
7258 }
7259
7260 double factorial (int start, int end);
7261 double M_chooses_N ( int m, int N)
7262 {
7263   /*Choose M elemets in N*/
7264   if ( m==N) return 1;
7265   else if ( m>N)
7266     {
7267       fprintf ( stderr, "\nERROR: M chosses N out of bounds ( M>N) [FATAL:%s]", PROGRAM);
7268       myexit (EXIT_FAILURE);
7269     }
7270   else if ( N<50)
7271     {
7272       return factorial (m+1, N)/factorial (1, N-m);
7273     }
7274   else
7275     {
7276       fprintf ( stderr, "\nERROR: M chosses N out of bounds ( N>50). Use log space [FATAL:%s]", PROGRAM);
7277       myexit (EXIT_FAILURE);
7278     }
7279   return -1;
7280 }
7281 double factorial (int start, int end)
7282   {
7283
7284     if ( start>end || start<0 || end<0)
7285       {
7286         fprintf ( stderr, "\nERROR: Negative Factorial [FATAL:%s]", PROGRAM);
7287         myexit ( EXIT_FAILURE);
7288       }
7289     else if (end==0) return 1;
7290     else if (end==start) return end;
7291     else
7292       {
7293         static double **lu;
7294         if ( !lu)lu=declare_double (100, 100);
7295
7296         if ( lu[start][end])return lu[start][end];
7297         else
7298           {
7299             int a;
7300             lu[start][end]=(double)start;
7301             for ( a=start+1; a<=end; a++)
7302               {
7303                 lu[start][end]*=(double)a;
7304               }
7305             return  lu[start][end];
7306           }
7307       }
7308     return -1;
7309   }
7310 /*********************************************************************/
7311 /*                                                                   */
7312 /*                         Fast Log Additions (adapted from Probcons)*/
7313 /*                                                                   */
7314 /*                                                                   */
7315 /*********************************************************************/
7316 double EXP (double x){
7317   //return exp(x);
7318   if (x > -2){
7319     if (x > -0.5){
7320       if (x > 0)
7321         return exp(x);
7322       return (((0.03254409303190190000*x + 0.16280432765779600000)*x + 0.49929760485974900000)*x + 0.99995149601363700000)*x + 0.99999925508501600000;
7323     }
7324     if (x > -1)
7325       return (((0.01973899026052090000*x + 0.13822379685007000000)*x + 0.48056651562365000000)*x + 0.99326940370383500000)*x + 0.99906756856399500000;
7326     return (((0.00940528203591384000*x + 0.09414963667859410000)*x + 0.40825793595877300000)*x + 0.93933625499130400000)*x + 0.98369508190545300000;
7327   }
7328   if (x > -8){
7329     if (x > -4)
7330       return (((0.00217245711583303000*x + 0.03484829428350620000)*x + 0.22118199801337800000)*x + 0.67049462206469500000)*x + 0.83556950223398500000;
7331     return (((0.00012398771025456900*x + 0.00349155785951272000)*x + 0.03727721426017900000)*x + 0.17974997741536900000)*x + 0.33249299994217400000;
7332   }
7333   if (x > -16)
7334     return (((0.00000051741713416603*x + 0.00002721456879608080)*x + 0.00053418601865636800)*x + 0.00464101989351936000)*x + 0.01507447981459420000;
7335   return 0;
7336 }
7337
7338 float LOOKUP (float x){
7339
7340   if (x <= 1.00f) return ((-0.009350833524763f * x + 0.130659527668286f) * x + 0.498799810682272f) * x + 0.693203116424741f;
7341   if (x <= 2.50f) return ((-0.014532321752540f * x + 0.139942324101744f) * x + 0.495635523139337f) * x + 0.692140569840976f;
7342   if (x <= 4.50f) return ((-0.004605031767994f * x + 0.063427417320019f) * x + 0.695956496475118f) * x + 0.514272634594009f;
7343
7344   return ((-0.000458661602210f * x + 0.009695946122598f) * x + 0.930734667215156f) * x + 0.168037164329057f;
7345 }
7346 void LOG_PLUS_EQUALS (float *x, float y){
7347
7348   if (x[0] < y)
7349     x[0] = (x[0] == LOG_ZERO || y - x[0] >= LOG_UNDERFLOW_THRESHOLD) ? y : LOOKUP(y-x[0]) + x[0];
7350   else
7351     x[0] = (y == LOG_ZERO || x[0] - y >= LOG_UNDERFLOW_THRESHOLD) ? x[0]  : LOOKUP(x[0]-y) + y;
7352 }
7353
7354 float LOG_ADD (float x, float y){
7355   if (x < y) return (x == LOG_ZERO || y - x >= LOG_UNDERFLOW_THRESHOLD) ? y : LOOKUP(y-x) + x;
7356   return (y == LOG_ZERO || x - y >= LOG_UNDERFLOW_THRESHOLD) ? x : LOOKUP(x-y) + y;
7357 }
7358
7359 float LOG_ADD3 (float x1, float x2, float x3){
7360   return LOG_ADD (x1, LOG_ADD (x2, x3));
7361 }
7362 float LOG_ADD4 (float x1, float x2, float x3, float x4){
7363   return LOG_ADD (x1, LOG_ADD (x2, LOG_ADD (x3, x4)));
7364 }
7365 float LOG_ADD5 (float x1, float x2, float x3, float x4, float x5){
7366   return LOG_ADD (x1, LOG_ADD (x2, LOG_ADD (x3, LOG_ADD (x4, x5))));
7367 }
7368 float LOG_ADD6 (float x1, float x2, float x3, float x4, float x5, float x6){
7369   return LOG_ADD (x1, LOG_ADD (x2, LOG_ADD (x3, LOG_ADD (x4, LOG_ADD (x5, x6)))));
7370 }
7371 float LOG_ADD7 (float x1, float x2, float x3, float x4, float x5, float x6, float x7){
7372   return LOG_ADD (x1, LOG_ADD (x2, LOG_ADD (x3, LOG_ADD (x4, LOG_ADD (x5, LOG_ADD (x6, x7))))));
7373 }
7374
7375
7376 #define LONG_SIZE 2
7377 #define SHORT_SIZE 1
7378 #define SPACE_PAD 4
7379 #define STD_SIZE 0
7380 char *strscn(char *s, char *pattern);
7381 long unsigned strtou(char *s, int base, char **scan_end);
7382 long int strtoi(char *s, int base, char **scan_end);
7383 int my_isnumber(char c, int base);
7384 int tonumber(char c);
7385
7386 int my_vsscanf(char *buf, char *fmt, va_list parms)
7387  {
7388          int scanned = 0, size = 0, suppress = 0;
7389          int w = 0, flag = 0, l = 0;
7390          char c, *c_ptr;
7391          long int n1, *n1l;
7392          int *n1b;
7393          short int *n1s;
7394          long unsigned n2, *n2l, parsing = 0;
7395          unsigned *n2b;
7396          short unsigned *n2s;
7397          double n3, *n3l;
7398          float *n3s;
7399          char *base = buf;
7400          while (*fmt != 0) {
7401                  if (*fmt != '%' && !parsing) {
7402                          /* No token detected */
7403                          fmt++;
7404                  } else {
7405                          /* We need to make a conversion */
7406                          if (*fmt == '%') {
7407                                  fmt++;
7408                                  parsing = 1;
7409                                  size = STD_SIZE;
7410                                  suppress = 0;
7411                                  w = 0;
7412                                  flag = 0;
7413                                  l = 0;
7414                          }
7415                          /* Parse token */
7416                          switch (*fmt) {
7417                          case '1':
7418                          case '2':
7419                          case '3':
7420                          case '4':
7421                          case '5':
7422                          case '6':
7423                          case '7':
7424                          case '8':
7425                          case '9':
7426                          case '0':
7427                                  if (parsing == 1) {
7428                                          w = strtou(fmt, 10, &base);
7429                                          /* We use SPACE_PAD to parse %10s
7430                                           * commands where the number is the
7431                                           * maximum number of char to store!
7432                                           */
7433                                          flag |= SPACE_PAD;
7434                                          fmt = base - 1;
7435                                  }
7436                                  break;
7437                          case 'c':
7438                                  c = *buf++;
7439                                  c_ptr = va_arg(parms, char *);
7440                                  *c_ptr = c;
7441                                  scanned++;
7442                                  parsing = 0;
7443                                  break;
7444                          case 's':
7445                                  c_ptr = va_arg(parms, char *);
7446                                  while (*buf != 0 && isspace(*buf))
7447                                          buf++;
7448                                  l = 0;
7449                                  while (*buf != 0 && !isspace(*buf)) {
7450                                          if (!(flag & SPACE_PAD))
7451                                                  *c_ptr++ = *buf;
7452                                          else if (l < w) {
7453                                                  *c_ptr++ = *buf;
7454                                                  l++;
7455                                          }
7456                                          buf++;
7457                                  }
7458                                  *c_ptr = 0;
7459                                  scanned++;
7460                                  parsing = 0;
7461                                  break;
7462                          case 'i':
7463                          case 'd':
7464                                  buf = strscn(buf, "1234567890-+");
7465                                  n1 = strtoi(buf, 10, &base);
7466                                  buf = base;
7467                                  if (!suppress) {
7468                                          switch (size) {
7469                                          case STD_SIZE:
7470                                                  n1b = va_arg(parms, int *);
7471                                                  *n1b = (int) n1;
7472                                                  break;
7473                                          case LONG_SIZE:
7474                                                  n1l = va_arg(parms,
7475                                                               long int *);
7476                                                  *n1l = n1;
7477                                                  break;
7478                                          case SHORT_SIZE:
7479                                                  n1s = va_arg(parms,
7480                                                               short int *);
7481                                                  *n1s = (short) (n1);
7482                                                  break;
7483                                          }
7484                                          scanned++;
7485                                  }
7486                                  parsing = 0;
7487                                  break;
7488                          case 'u':
7489                                  buf = strscn(buf, "1234567890");
7490                                  n2 = strtou(buf, 10, &base);
7491                                  buf = base;
7492                                  if (!suppress) {
7493                                          switch (size) {
7494                                          case STD_SIZE:
7495                                                  n2b = va_arg(parms,
7496                                                               unsigned *);
7497                                                  *n2b = (unsigned) n2;
7498                                                  break;
7499                                          case LONG_SIZE:
7500                                                  n2l = va_arg(parms,
7501                                                               long unsigned *);
7502                                                  *n2l = n2;
7503                                                  break;
7504                                          case SHORT_SIZE:
7505                                                  n2s = va_arg(parms, short unsigned
7506                                                               *);
7507                                                  *n2s = (short) (n2);
7508                                                  break;
7509                                          }
7510                                          scanned++;
7511                                  }
7512                                  parsing = 0;
7513                                  break;
7514                          case 'x':
7515                                  buf = strscn(buf, "1234567890xabcdefABCDEF");
7516                                  n2 = strtou(buf, 16, &base);
7517                                  buf = base;
7518                                  if (!suppress) {
7519                                          switch (size) {
7520                                          case STD_SIZE:
7521                                                  n2b = va_arg(parms,
7522                                                               unsigned *);
7523                                                  *n2b = (unsigned) n2;
7524                                                  break;
7525                                          case LONG_SIZE:
7526                                                  n2l = va_arg(parms,
7527                                                               long unsigned *);
7528                                                  *n2l = n2;
7529                                                  break;
7530                                          case SHORT_SIZE:
7531                                                  n2s = va_arg(parms, short unsigned
7532                                                               *);
7533                                                  *n2s = (short) (n2);
7534                                                  break;
7535                                          }
7536                                          scanned++;
7537                                  }
7538                                  parsing = 0;
7539                                  break;
7540                          case 'f':
7541                          case 'g':
7542                          case 'e':
7543                                  buf = strscn(buf, "1234567890.e+-");
7544                                  n3 = strtod(buf, &base);
7545                                  buf = base;
7546                                  if (!suppress) {
7547                                          switch (size) {
7548                                          case STD_SIZE:
7549                                                  n3l = va_arg(parms, double *);
7550                                                  *n3l = n3;
7551                                                  break;
7552                                          case LONG_SIZE:
7553                                                  n3l = va_arg(parms, double *);
7554                                                  *n3l = n3;
7555                                                  break;
7556                                          case SHORT_SIZE:
7557                                                  n3s = va_arg(parms, float *);
7558                                                  *n3s = (float) (n3);
7559                                                  break;
7560                                          }
7561                                          scanned++;
7562                                  }
7563                                  parsing = 0;
7564                                  break;
7565                          case 'l':
7566                                  size = LONG_SIZE;
7567                                  break;
7568                          case 'h':
7569                          case 'n':
7570                                  size = SHORT_SIZE;
7571                                  break;
7572                          case '*':
7573                                  suppress = 1;
7574                                  break;
7575                          default:
7576                                  parsing = 0;
7577                                  break;
7578                          }
7579                          fmt++;
7580                  }
7581          }
7582          return (scanned);
7583  }
7584 char *strscn(char *s, char *pattern)
7585  {
7586          char *scan;
7587          while (*s != 0) {
7588                  scan = pattern;
7589                  while (*scan != 0) {
7590                          if (*s == *scan)
7591                                  return (s);
7592                          else
7593                                  scan++;
7594                  }
7595                  s++;
7596          }
7597          return (NULL);
7598  }
7599
7600 long unsigned strtou(char *s, int base, char **scan_end)
7601  {
7602          int value, overflow = 0;
7603          long unsigned result = 0, oldresult;
7604          /* Skip trailing zeros */
7605          while (*s == '0')
7606                  s++;
7607          if (*s == 'x' && base == 16) {
7608                  s++;
7609                  while (*s == '0')
7610                          s++;
7611          }
7612          /* Convert number */
7613          while (my_isnumber(*s, base)) {
7614                  value = tonumber(*s++);
7615                  if (value > base || value < 0)
7616                          return (0);
7617                  oldresult = result;
7618                  result *= base;
7619                  result += value;
7620                  /* Detect overflow */
7621                  if (oldresult > result)
7622                          overflow = 1;
7623          }
7624          if (scan_end != 0L)
7625                  *scan_end = s;
7626          if (overflow)
7627                  result = INT_MAX;
7628          return (result);
7629  }
7630 long int strtoi(char *s, int base, char **scan_end)
7631  {
7632          int sign, value, overflow = 0;
7633          long int result = 0, oldresult;
7634          /* Evaluate sign */
7635          if (*s == '-') {
7636                  sign = -1;
7637                  s++;
7638          } else if (*s == '+') {
7639                  sign = 1;
7640                  s++;
7641          } else
7642                  sign = 1;
7643          /* Skip trailing zeros */
7644          while (*s == '0')
7645                  s++;
7646          /* Convert number */
7647          while (my_isnumber(*s, base)) {
7648                  value = tonumber(*s++);
7649                  if (value > base || value < 0)
7650                          return (0);
7651                  oldresult = result;
7652                  result *= base;
7653                  result += value;
7654                  /* Detect overflow */
7655                  if (oldresult > result)
7656                          overflow = 1;
7657          }
7658          if (scan_end != 0L)
7659                  *scan_end = s;
7660          if (overflow)
7661                  result = INT_MAX;
7662          result *= sign;
7663          return (result);
7664  }
7665
7666 int my_isnumber(char c, int base)
7667  {
7668          static char *digits = "0123456789ABCDEF";
7669          if ((c >= '0' && c <= digits[base - 1]))
7670                  return (1);
7671          else
7672                  return (0);
7673  }
7674
7675  int tonumber(char c)
7676  {
7677          if (c >= '0' && c <= '9')
7678                  return (c - '0');
7679          else if (c >= 'A' && c <= 'F')
7680                  return (c - 'A' + 10);
7681          else if (c >= 'a' && c <= 'f')
7682                  return (c - 'a' + 10);
7683          else
7684                  return (c);
7685  }
7686
7687 ///////////////////////////////////////////////////////////////////////////////////////////
7688 // Hash function
7689 ////////////////////////////////////////////////////////////////////////////////////////////
7690 unsigned long hash_file(char* file)  //returns the hash value for key
7691     {
7692       // Calculate a hash value by the division method:
7693       // Transform key into a natural number k = sum ( key[i]*128^(L-i) ) and calculate i= k % num_slots.
7694       // Since calculating k would lead to an overflow, i is calculated iteratively
7695       // and at each iteration the part divisible by num_slots is subtracted, i.e. (% num_slots is taken).
7696
7697       unsigned long i=0;     // Start of iteration: k is zero
7698       unsigned long num_slots=999999999;
7699
7700
7701       FILE *fp;
7702       unsigned long c;
7703
7704
7705       if (file==NULL || !check_file_exists (file) ) {printf("Warning from util.c:hasch_file: No File [FATAL:%s]\n", PROGRAM); myexit (EXIT_FAILURE);}
7706       num_slots/=128;
7707       fp=vfopen (file, "r");
7708       while ( (c=fgetc (fp))!=EOF)
7709         {
7710           i = ((i<<7) + c) % num_slots;
7711         }
7712       vfclose (fp);
7713
7714       return i;
7715     }
7716 int ** r_generate_array_int_list ( int len, int min, int max,int step, int **array, int f, int *n,FILE *fp, int *c_array);
7717 int **generate_array_int_list (int len, int min, int max, int step, int *n, char *file)
7718    {
7719      int **array, *c_array;
7720      FILE *fp=NULL;
7721
7722      if (n==NULL)
7723        {
7724          array=NULL;
7725          fp=vfopen (file, "w");
7726        }
7727      else
7728        {
7729          int a,s;
7730          n[0]=0;
7731          for (s=1, a=0; a<len; a++)s*=((max-min)+1)/step;
7732          array=declare_int (s, len+1);
7733        }
7734      c_array=vcalloc (len, sizeof (int));
7735      array=r_generate_array_int_list ( len, max, min, step, array, 0, n,fp, c_array);
7736      vfree (c_array);
7737      if ( fp) vfclose (fp);
7738      return array;
7739    }
7740 int ** r_generate_array_int_list ( int len, int min, int max,int step, int **array, int f, int *n,FILE *fp, int *c_array)
7741    {
7742      int a;
7743
7744      if ( f==len)
7745        {
7746          if ( array)
7747            {
7748              for (a=0; a<len; a++)
7749                {
7750                  array[n[0]][a]=c_array[a];
7751                }
7752              n[0]++;
7753            }
7754          else
7755            {
7756
7757
7758              for (a=0; a<len; a++)fprintf ( fp, "%3d ",c_array[a]);
7759              fprintf (fp, "\n");
7760            }
7761          return array;
7762
7763        }
7764      else
7765        {
7766          for (a=max; a<=min;a+=step)
7767            {
7768              c_array[f]=a;
7769              r_generate_array_int_list (len, min, max, step, array, f+1, n,fp, c_array);
7770            }
7771        }
7772      return array;
7773    }
7774
7775 char *** r_generate_array_string_list ( int len, char ***alp,int *alp_size, char ***array, int f, int *n,FILE *fp, char **c_array, int mode, int pstart);
7776 char ***generate_array_string_list (int len, char ***alp, int *alp_size, int *n, char *file, int mode)
7777    {
7778      char  ***array, **c_array;
7779      FILE *fp=NULL;
7780
7781      if (file!=NULL)
7782        {
7783          array=NULL;
7784          n[0]=0;
7785          fp=vfopen (file, "w");
7786        }
7787      else
7788        {
7789
7790          int a,s;
7791          n[0]=0;
7792          for (s=1, a=0; a<len; a++)
7793            {
7794              s*=alp_size[a];
7795            }
7796          array=declare_arrayN (3,sizeof (char), s, len,0);
7797        }
7798      c_array=declare_char (len,0);
7799      array=r_generate_array_string_list ( len, alp, alp_size, array, 0, n,fp, c_array, mode, -1);
7800      vfree (c_array);
7801      if ( fp) vfclose (fp);
7802      return array;
7803    }
7804 char *** r_generate_array_string_list ( int len, char ***alp, int *alp_size, char  ***array, int f, int *n,FILE *fp, char **c_array, int mode, int pstart)
7805    {
7806      int a;
7807      int start;
7808
7809      if ( f==len)
7810        {
7811          if ( array)
7812            {
7813              for (a=0; a<len; a++)
7814                {
7815                  array[n[0]][a]=c_array[a];
7816                }
7817              n[0]++;
7818            }
7819          else
7820            {
7821
7822              n[0]++;
7823              for (a=0; a<len; a++)fprintf ( fp, "%s ",c_array[a]);
7824              fprintf (fp, "\n");
7825            }
7826          return array;
7827
7828        }
7829      else
7830        {
7831          if ( mode==OVERLAP)
7832            {
7833              start=0;
7834            }
7835          else if ( mode==NO_OVERLAP)
7836            {
7837              start=pstart+1;
7838            }
7839
7840          for (a=start; a<alp_size[f]; a++)
7841            {
7842              c_array[f]=alp[f][a];
7843              r_generate_array_string_list (len,alp, alp_size, array, f+1, n,fp, c_array, mode, a);
7844            }
7845        }
7846      return array;
7847    }
7848
7849 float rates2sensitivity (int tp, int tn, int fp, int fn, float *sp, float *sn, float *sen2, float *b)
7850 {
7851   if (sp==NULL)
7852     {
7853       sp=vcalloc (1, sizeof (float));
7854       sn=vcalloc (1, sizeof (float));
7855       sen2=vcalloc (1, sizeof (float));
7856       b=vcalloc (1, sizeof (float));
7857     }
7858   sn[0]  =((tp+fn)==0)?1:(float)tp/(float)(tp+fn);
7859   sen2[0]=((tp+fp)==0)?1:(float)tp/(float)(tp+fp);
7860   sp[0]  =((tn+fp)==0)?1:(float)tn/(float)(tn+fp);
7861   b[0]=MIN((MIN((sp[0]),(sn[0]))),(sen2[0]));
7862   return b[0];
7863 }
7864 float profile2sensitivity (char *pred, char *ref, float *sp, float *sn, float *sen2, float *b)
7865 {
7866   int tp=0, tn=0, fp=0, fn=0;
7867   int a, l;
7868
7869   l=strlen (pred);
7870
7871   for (a=0; a<l; a++)
7872     {
7873       if (pred[a]=='I' && ref[a]=='I')tp++;
7874       if (pred[a]=='O' && ref[a]=='I')fn++;
7875       if (pred[a]=='I' && ref[a]=='O')fp++;
7876       if (pred[a]=='O' && ref[a]=='O')tn++;
7877     }
7878   return  rates2sensitivity (tp, tn, fp, fn, sp, sn, sen2, b);
7879 }
7880 float profile2evalue (char *pred, char *ref)
7881 {
7882
7883   int a, l;
7884   double P=0;
7885   double E=0;
7886   double II=0;
7887   double p1, p2, p;
7888   l=strlen (pred);
7889
7890   for (a=0; a<l; a++)
7891     {
7892       if (pred[a]=='I')P++;
7893       if (ref[a]=='I') E++;
7894
7895       if (pred[a]=='I' && ref[a]=='I')II++;
7896     }
7897
7898   if (II==0)return 0;
7899
7900   p1= M_chooses_Nlog (P,l) + M_chooses_Nlog (II, P) + M_chooses_Nlog (E-II, l-P);
7901   p2=(M_chooses_Nlog (P,l)+M_chooses_Nlog (E,l));
7902   p=(p1-p2);
7903
7904   return (float)-p;
7905 }
7906
7907
7908
7909
7910
7911 /*********************************COPYRIGHT NOTICE**********************************/
7912 /*� Centro de Regulacio Genomica */
7913 /*and */
7914 /*Cedric Notredame */
7915 /*Sat Aug 16 20:30:21 WEST 2008. */
7916 /*All rights reserved.*/
7917 /*This file is part of T-COFFEE.*/
7918 /**/
7919 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
7920 /*    it under the terms of the GNU General Public License as published by*/
7921 /*    the Free Software Foundation; either version 2 of the License, or*/
7922 /*    (at your option) any later version.*/
7923 /**/
7924 /*    T-COFFEE is distributed in the hope that it will be useful,*/
7925 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
7926 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
7927 /*    GNU General Public License for more details.*/
7928 /**/
7929 /*    You should have received a copy of the GNU General Public License*/
7930 /*    along with Foobar; if not, write to the Free Software*/
7931 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
7932 /*...............................................                                                                                      |*/
7933 /*  If you need some more information*/
7934 /*  cedric.notredame@europe.com*/
7935 /*...............................................                                                                                                                                     |*/
7936 /**/
7937 /**/
7938 /*      */
7939 /*********************************COPYRIGHT NOTICE**********************************/
7940 /*
7941 FILE * farray()
7942 {
7943   return tmpfile();
7944 }
7945 float fset (FILE *fp, int i, float v)
7946 {
7947   fseek (fp, i*sizeof (float), SEEK_SET);
7948   fwrite(&v, sizeof (float), fp);
7949   return v;
7950 }
7951 float fread (FILE *fp, int i, float v)
7952 {
7953   fseek (fp, i*sizeof (float), SEEK_SET);
7954   fread (&v, sizeof (float), fp);
7955   return v;
7956 }
7957 */
7958 /*********************************COPYRIGHT NOTICE**********************************/
7959 /*© Centro de Regulacio Genomica */
7960 /*and */
7961 /*Cedric Notredame */
7962 /*Tue Oct 27 10:12:26 WEST 2009. */
7963 /*All rights reserved.*/
7964 /*This file is part of T-COFFEE.*/
7965 /**/
7966 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
7967 /*    it under the terms of the GNU General Public License as published by*/
7968 /*    the Free Software Foundation; either version 2 of the License, or*/
7969 /*    (at your option) any later version.*/
7970 /**/
7971 /*    T-COFFEE is distributed in the hope that it will be useful,*/
7972 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
7973 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
7974 /*    GNU General Public License for more details.*/
7975 /**/
7976 /*    You should have received a copy of the GNU General Public License*/
7977 /*    along with Foobar; if not, write to the Free Software*/
7978 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
7979 /*...............................................                                                                                      |*/
7980 /*  If you need some more information*/
7981 /*  cedric.notredame@europe.com*/
7982 /*...............................................                                                                                                                                     |*/
7983 /**/
7984 /**/
7985 /*      */
7986 /*********************************COPYRIGHT NOTICE**********************************/