binaries/src/clustalo/src/hhalign/util-C.h

   1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
   2
   3 /*********************************************************************
   4  * Clustal Omega - Multiple sequence alignment
   5  *
   6  * Copyright (C) 2010 University College Dublin
   7  *
   8  * Clustal-Omega is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License as
  10  * published by the Free Software Foundation; either version 2 of the
  11  * License, or (at your option) any later version.
  12  *
  13  * This file is part of Clustal-Omega.
  14  *
  15  ********************************************************************/
  16
  17 /*
  18  * RCS $Id: util-C.h 155 2010-11-17 12:18:47Z fabian $
  19  */
  20
  21 // Utility subroutines
  22
  23
  24 #ifndef MAIN
  25 #include <iostream>   // cin, cout, cerr
  26 #include <fstream>    // ofstream, ifstream
  27 #include <cstdio>     // printf
  28 #include <stdlib.h>   // exit
  29 #include <time.h>     // clock
  30 #endif
  31 #include <sys/time.h>
  32
  33 /////////////////////////////////////////////////////////////////////////////////////
  34 // Arithmetics
  35 /////////////////////////////////////////////////////////////////////////////////////
  36
  37 //// max and min
  38 inline double dmax(double x, double y) { return (x>y? x : y);}
  39 inline double dmin(double x, double y) { return (x<y? x : y);}
  40 inline int imax(int x, int y) { return (x>y? x : y);}
  41 inline int imin(int x, int y) { return (x<y? x : y);}
  42 inline int iabs(int x) { return (x>=0? x : -x);}
  43
  44 // Rounding up, rounding down and rounding to nearest integer
  45 inline int iceil(double x)  {return int(ceil(x));}
  46 inline int ifloor(double x) {return int(floor(x));}
  47 inline int iround(double x) {return int(floor(x+0.5));}
  48
  49 //// Generalized mean: d=0: sqrt(x*y)  d=1: (x+y)/2  d->-inf: min(x,y)  d->+inf: max(x,y)
  50 inline double fmean(double x, double y, double d) { return pow( (pow(x,d)+pow(y,d))/2 ,1./d);}
  51
  52 // log base 2
  53 inline float log2(float x)  {return (x<=0? (float)(-100000):1.442695041*log(x));}
  54 inline float log10(float x) {return (x<=0? (float)(-100000):0.434294481*log(x));}
  55
  56
  57 /////////////////////////////////////////////////////////////////////////////////////
  58 // fast log base 2
  59 /////////////////////////////////////////////////////////////////////////////////////
  60
  61 // This function returns log2 with a max abolute deviation of +/- 1.5E-5 (typically 0.8E-5).
  62 // It takes 1.42E-8 s  whereas log2(x) takes 9.5E-7 s. It is hence 9.4 times faster.
  63 // It makes use of the representation of 4-byte floating point numbers:
  64 // seee eeee emmm mmmm mmmm mmmm mmmm mmmm
  65 // s is the sign,
  66 // the following 8 bits, eee eee e, give the exponent + 127 (in hex: 0x7f).
  67 // The following 23 bits, m, give the mantisse, the binary digits behind the decimal point.
  68 // In summary: x = (-1)^s * 1.mmmmmmmmmmmmmmmmmmmmmm * 2^(eeeeeee-127)
  69 // The expression (((*(int *)&x) & 0x7f800000 ) >>23 )-0x7f is the exponent eeeeeeee, i.e.
  70 // the largest integer that is smaller than log2(x) (e.g. -1 for 0.9). *(int *)&x is an integer which
  71 // contains the bytes as the floating point variable x is represented in memory.
  72 // Check:  assert( sizeof(f) == sizeof(int) );
  73 // Check:  assert( sizeof(f) == 4 );
  74 inline float fast_log2(float x)
  75 {
  76   static float lg2[1025];         // lg2[i] = log2[1+x/1024]
  77   static float diff[1025];        // diff[i]= (lg2[i+1]-lg2[i])/8096 (for interpolation)
  78   static char initialized;
  79   if (x<=0) return -100000;
  80   if (!initialized)   //First fill in the arrays lg2[i] and diff[i]
  81     {
  82       float prev = 0.0f;
  83       lg2[0] = 0.0f;
  84       for (int i=1; i<=1024; ++i)
  85         {
  86           lg2[i] = log(float(1024+i))*1.442695041-10.0f;
  87           diff[i-1] = (lg2[i]-prev)*1.2352E-4;
  88           prev = lg2[i];
  89         }
  90       initialized=1;
  91     }
  92   int a = (((*((int *)&x)) & 0x7F800000) >>23 )-0x7f;
  93   int b =  ((*((int *)&x)) & 0x007FE000) >>13;
  94   int c =  ((*((int *)&x)) & 0x00001FFF);
  95   return a + lg2[b] + diff[b]*(float)(c);
  96 }
  97
  98 /////////////////////////////////////////////////////////////////////////////////////
  99 // fast 2^x
 100 // ATTENTION: need to compile with g++ -fno-strict-aliasing when using -O2 or -O3!!!
 101 // Relative deviation < 1.5E-4
 102 /////////////////////////////////////////////////////////////////////////////////////
 103 inline float fpow2(float x)
 104 {
 105   if (x>=128) return FLT_MAX;
 106   if (x<=-128) return FLT_MIN;
 107   int *px = (int*)(&x);                 // store address of float as pointer to long
 108   float tx = (x-0.5f) + (3<<22);        // temporary value for truncation: x-0.5 is added to a large integer (3<<22)
 109   int lx = *((int*)&tx) - 0x4b400000;   // integer value of x
 110   float dx = x-(float)(lx);             // float remainder of x
 111   x = 1.0f + dx*(0.6960656421638072f          // cubic apporoximation of 2^x
 112            + dx*(0.224494337302845f           // for x in the range [0, 1]
 113            + dx*(0.07944023841053369f)));
 114   *px += (lx<<23);                            // add integer power of 2 to exponent
 115   return x;
 116 }
 117
 118 /////////////////////////////////////////////////////////////////////////////////////
 119 // ATTENTION:
 120 // Can't be used with -O2/-O3 optimization on some compilers !
 121 // Works with g++ version 4.1, but not with 3.4, in which case it returns values
 122 // that are a factor 1.002179942 too low
 123 //
 124 // Fast pow2 routine (Johannes Soeding)
 125 // Same speed as fpow2(), but *relative* deviation < 1.2E-7
 126 // Makes use of the binary representation of floats in memory:
 127 //   x = (-1)^s * 1.mmmmmmmmmmmmmmmmmmmmmm * 2^(eeeeeee-127)
 128 // is represented as
 129 // 31        23                   7654 3210
 130 //  seee eeee emmm mmmm mmmm mmmm mmmm mmmm
 131 // s is the sign, the 8 bits eee eee e are the exponent + 127 (in hex: 0x7f),
 132 // and the following 23 bits m give the mantisse.
 133 // We decompose the argument x = a + b, with integer a and 0 <= b < 1
 134 // Therefore 2^x = 2^a * 2^b  where a is the binary exponent of 2^x
 135 // and  1 <= 2^b < 2,  i.e. 2^b determines the mantisse uniquely.
 136 // To calculate 2^b, we split b into the first 10 bits and the last 13 bits,
 137 // b = b' + c, and then look up the mantisse of 2^b' in a precomputed table.
 138 // We use the residual c to interpolate between the mantisse for 2^b' and 2(b'+1/1024)
 139 /////////////////////////////////////////////////////////////////////////////////////
 140 inline float fast_pow2(float x)
 141 {
 142   if (x<=-127) return 5.9E-39;
 143   if (x>=128)  return 3.4E38;
 144   static char initialized=0;
 145   static unsigned int pow2[1025];
 146   static unsigned int diff[1025];
 147   static int y = 0;
 148   if (!initialized)   //First fill in the pow2-vector
 149     {
 150       float f;
 151       unsigned int prev = 0;
 152       pow2[0] = 0;
 153       for (int b=1; b<1024; b++)
 154         {
 155           f=pow(2.0,float(b)/1024.0);
 156           pow2[b]=(*((unsigned int *)(&f)) & 0x7FFFFF); // store the mantisse of 2^(1+b/1024)
 157           diff[b-1]=pow2[b]-prev;
 158           prev=pow2[b];
 159         }
 160       pow2[1024]=0x7FFFFF;
 161       diff[1023]=pow2[1024]-prev;
 162       initialized=1;
 163     }
 164
 165   int *px = (int *)(&x);                              // store address of float as pointer to int
 166   int E = ((*px & 0x7F800000)>>23)-127;               // E is exponent of x and is <=6
 167   unsigned int M=(*px & 0x007FFFFF) | 0x00800000;     // M is the mantisse 1.mmm mmmm mmmm mmmm mmmm mmmm
 168   int a,b,c;
 169   if (x>=0)
 170     {
 171       if (E>=0) {
 172         a = 0x3F800000 + ((M<<E) & 0x7F800000);       // a is exponent of 2^x, beginning at bit 23
 173         b = ((M<<E) & 0x007FE000)>>13;
 174         c = ((M<<E) & 0x00001FFF);
 175       } else {
 176         a =  0x3F800000;                              // a = exponent of 2^x = 0
 177         b = ((M>>(-E)) & 0x007FE000)>>13;
 178         c = ((M>>(-E)) & 0x00001FFF);
 179       }
 180     }
 181   else
 182     {
 183       if (E>=0) {
 184         a = 0x3F000000 - ((M<<E) & 0x7F800000);       // a is exponent of 2^x
 185         b = (0x00800000-(int)((M<<E) & 0x007FFFFF)) >>13;
 186         c = (0x00800000-(int)((M<<E) & 0x007FFFFF)) & 0x00001FFF;
 187       } else {
 188         a = 0x3F000000;                               // a = exponent of 2^x = -1
 189         b = (0x00800000-(int)((M>>(-E)) & 0x007FFFFF)) >>13;
 190         c = (0x00800000-(int)((M>>(-E)) & 0x007FFFFF)) & 0x00001FFF;
 191       }
 192     }
 193 /*   printf("x=%0X\n",*px); */
 194 /*   printf("E=%0X\n",E); */
 195 /*   printf("M=%0X\n",M); */
 196 /*   printf("a=%0X\n",a); */
 197 /*   printf("b=%0X\n",b); */
 198   y = a | (pow2[b] + ((diff[b]*c)>>13) );
 199   /*   printf("2^x=%0X\n",*px); */
 200   return *((float*)&y);
 201 }
 202
 203
 204
 205 // Normalize a float array such that it sums to one
 206 // If it sums to 0 then assign def_array elements to array (optional)
 207 inline float NormalizeTo1(float* array, int length, float* def_array=NULL)
 208 {
 209   float sum=0.0f;
 210   int k;
 211   for (k=0; k<length; k++) sum+=array[k];
 212   if (sum!=0.0f)
 213     {
 214       float fac=1.0/sum;
 215       for (k=0; k<length; k++) array[k]*=fac;
 216     }
 217   else if (def_array)
 218     for (k=0; k<length; k++) array[k]=def_array[k];
 219   return sum;
 220 }
 221
 222 // Normalize a float array such that it sums to x
 223 // If it sums to 0 then assign def_array elements to array (optional)
 224 inline float NormalizeToX(float* array, int length, float x, float* def_array=NULL)
 225 {
 226   float sum=0.0;
 227   int k;
 228   for (k=0; k<length; k++) sum+=array[k];
 229   if (sum)
 230     {
 231       float fac=x/sum;
 232       for (k=0; k<length; k++) array[k]*=fac;
 233     }
 234   else if (def_array)
 235     for (k=0; k<length; k++) array[k]=def_array[k];
 236   return sum;
 237 }
 238
 239 /////////////////////////////////////////////////////////////////////////////////////
 240 // Similar to spintf("%*g,w,val), but displays maximum number of digits within width w
 241 /////////////////////////////////////////////////////////////////////////////////////
 242 inline char* sprintg(float val, int w)
 243 {
 244   static char str[100];
 245   float log10val = log10(fabs(val));
 246   int neg = (val<0? 1: 0);
 247   if (log10val >= w-neg-1 || -log10val > 3)
 248     {
 249       // positive exponential 1.234E+06
 250       // negative exponential 1.234E-06
 251       int d = w-6-neg;
 252       sprintf(str,"%*.*e",w,d<1?1:d,val);
 253     }
 254   else
 255     {
 256       int d = log10val>0? w-2-neg-int(log10val): w-2-neg;
 257       sprintf(str,"%#*.*f",w,d,val);
 258     }
 259   return str;
 260 }
 261
 262 /////////////////////////////////////////////////////////////////////////////////////
 263 // String utilities
 264 /////////////////////////////////////////////////////////////////////////////////////
 265
 266 //the integer. If no integer is found, returns INT_MIN and sets ptr to NULL      /* MR1 */
 267 inline int strtoi(const char*& ptr)
 268 {
 269     int i;
 270     const char* ptr0=ptr;
 271     if (!ptr) return INT_MIN;
 272     while (*ptr!='\0' && !(*ptr>='0' && *ptr<='9')) ptr++;
 273     if (*ptr=='\0') {
 274         ptr=0;
 275         return INT_MIN;
 276     }
 277     if (*(ptr-1)=='-' && ptr>ptr0) i=-atoi(ptr); else i=atoi(ptr);
 278     while (*ptr>='0' && *ptr<='9') ptr++;
 279     return i;
 280 }
 281
 282
 283 //Same as strint, but interpretes '*' as default /* MR1 */
 284 inline int strtoi_(const char*& ptr, int deflt=INT_MAX)
 285 {
 286     int i;
 287     if (!ptr) return INT_MIN;
 288     while (*ptr!='\0' && !(*ptr>='0' && *ptr<='9') && *ptr!='*') ptr++;
 289     if (*ptr=='\0') {
 290         ptr=0;
 291         return INT_MIN;
 292     }
 293     if (*ptr=='*') {
 294         ptr++;
 295         return deflt;
 296     }
 297     if (*(ptr-1)=='-') i=atoi(ptr-1);
 298     else i=atoi(ptr);
 299     while (*ptr>='0' && *ptr<='9') ptr++;
 300     return i;
 301 }
 302
 303
 304 // Returns leftmost integer in ptr and sets the pointer to first char after
 305 // the integer. If no integer is found, returns INT_MIN and sets pt to NULL
 306 int strint(char*& ptr)
 307 {
 308   int i;
 309   char* ptr0=ptr;
 310   if (!ptr) return INT_MIN;
 311   while (*ptr!='\0' && !(*ptr>='0' && *ptr<='9')) ptr++;
 312   if (*ptr=='\0')
 313     {
 314       ptr=0;
 315       return INT_MIN;
 316     }
 317   if (*(ptr-1)=='-' && ptr>ptr0) i=-atoi(ptr); else i=atoi(ptr);
 318   while (*ptr>='0' && *ptr<='9') ptr++;
 319   return i;
 320 }
 321
 322 // Same as strint, but interpretes '*' as default
 323 int strinta(char*& ptr, int deflt=99999)
 324 {
 325   int i;
 326   if (!ptr) return INT_MIN;
 327   while (*ptr!='\0' && !(*ptr>='0' && *ptr<='9') && *ptr!='*') ptr++;
 328   if (*ptr=='\0')
 329     {
 330       ptr=0;
 331       return INT_MIN;
 332     }
 333   if (*ptr=='*')
 334     {
 335       ptr++;
 336       return deflt;
 337     }
 338   if (*(ptr-1)=='-') i=atoi(ptr-1);
 339   else i=atoi(ptr);
 340   while (*ptr>='0' && *ptr<='9') ptr++;
 341   return i;
 342 }
 343
 344 // Returns leftmost float in ptr and sets the pointer to first char after
 345 // the float. If no float is found, returns FLT_MIN and sets pt to NULL /* MR1 */
 346 float strflt(char*& ptr)
 347 {
 348     float i;
 349     char* ptr0=ptr;
 350     if (!ptr) return FLT_MIN;
 351     while (*ptr!='\0' && !(*ptr>='0' && *ptr<='9')) ptr++;
 352     if (*ptr=='\0')
 353         {
 354             ptr=0;
 355             return FLT_MIN;
 356         }
 357     if (ptr>ptr0 && *(ptr-1)=='-') i=-atof(ptr); else i=atof(ptr);
 358     while ((*ptr>='0' && *ptr<='9') || *ptr=='.') ptr++;
 359     return i;
 360 }
 361
 362 // Same as strint, but interpretes '*' as default  /* MR1 */
 363 float strflta(char*& ptr, float deflt=99999)
 364 {
 365     float i;
 366     if (!ptr) return FLT_MIN;
 367     while (*ptr!='\0' && !(*ptr>='0' && *ptr<='9') && *ptr!='*') ptr++;
 368     if (*ptr=='\0')
 369         {
 370             ptr=0;
 371             return FLT_MIN;
 372         }
 373     if (*ptr=='*')
 374         {
 375             ptr++;
 376             return deflt;
 377         }
 378     if (*(ptr-1)=='-') i=-atof(ptr);
 379     else i=atof(ptr);
 380     while ((*ptr>='0' && *ptr<='9') || *ptr=='.') ptr++;
 381     return i;
 382 }
 383
 384
 385 // Removes the newline and other control characters at the end of a string (if present)
 386 // and returns the new length of the string (-1 if str is NULL)
 387 inline int chomp(char str[])
 388 {
 389   if (!str) return -1;
 390   int l=0;
 391   for (l=strlen(str)-1; l>=0 && str[l]<32; l--);
 392   str[++l]='\0';
 393   return l;
 394 }
 395
 396 // Emulates the ifstream::getline method; similar to fgets(str,maxlen,FILE*),
 397 // but removes the newline at the end and returns NULL if at end of file or read error
 398 inline char* fgetline(char str[], const int maxlen, FILE* file)
 399 {
 400   if (!fgets(str,maxlen,file)) return NULL;
 401   if (chomp(str)+1>=maxlen)    // if line is cut after maxlen characters...
 402     while (fgetc(file)!='\n'); // ... read in rest of line
 403   return(str);
 404 }
 405
 406 // copies substring str[a,b] into substr and returns substr
 407 char *substr(char* substr, char* str, int a, int b)
 408 {
 409   if (b<a) {int i=b; b=a; a=i;}
 410   if (b-a>1000)
 411     {printf("Function substr: >1000 chars to copy. Exiting.\n"); exit(6);}
 412   char* dest=substr;
 413   char* source=str+a;
 414   char* send=str+b;
 415   while (*source!='\0' && source<=send) *(dest++) = *(source++);
 416   *dest='\0';
 417   return substr;
 418 }
 419
 420
 421 // Returns pointer to first non-white-space character in str OR to NULL if none found
 422 inline char* strscn(char* str)
 423 {
 424   if (!str) return NULL;
 425   char* ptr=str;
 426   while (*ptr!='\0' && *ptr<=32) ptr++;
 427   return (*ptr=='\0')? NULL: ptr;
 428 }
 429
 430 // Returns pointer to first white-space character in str OR to NULL if none found   /* MR1 */
 431 inline char* strscn_ws(char* str)
 432 {
 433     if (!str) return NULL;
 434     char* ptr=str;
 435     while (*ptr!='\0' && *ptr>32) ptr++;
 436     return (*ptr=='\0')? NULL: ptr;
 437 }
 438
 439 //Returns pointer to first non-white-space character in str OR to NULL if none found  /* MR1 */
 440 inline const char* strscn_c(const char* str)
 441 {
 442     if (!str) return NULL;
 443     const char* ptr=str;
 444     while (*ptr!='\0' && isspace(*ptr)) ptr++;
 445     return (*ptr=='\0') ? NULL : ptr;
 446 }
 447
 448 // Returns pointer to first  non-white-space character in str OR to end of string '\0' if none found
 449 inline char* strscn_(char* str)
 450 {
 451   if (!str) return NULL;
 452   char* ptr=str;
 453   while (*ptr!='\0' && *ptr<=32) ptr++;
 454   return ptr;
 455 }
 456
 457 // Returns pointer to first non-c character in str OR to NULL if none found
 458 inline char* strscn(char* str, const char c)
 459 {
 460   if (!str) return NULL;
 461   char* ptr=str;
 462   while (*ptr!='\0' && *ptr==c) ptr++;
 463   return (*ptr=='\0')? NULL: ptr;
 464 }
 465
 466 // Returns pointer to first  non-c character in str OR to end of string '\0' if none found
 467 inline char* strscn_(char* str, const char c)
 468 {
 469   if (!str) return NULL;
 470   char* ptr=str;
 471   while (*ptr!='\0' && *ptr==c) ptr++;
 472   return ptr;
 473 }
 474
 475 // Cuts string at first white space character found by overwriting it with '\0'.
 476 // Returns pointer to next non-white-space char OR to NULL if no such char found
 477 inline char* strcut(char* str)
 478 {
 479   if (!str) return NULL;
 480   char* ptr=str;
 481   while (*ptr!='\0' && *ptr>32) ptr++;
 482   if (*ptr=='\0') return NULL;
 483   *ptr='\0';
 484   ptr++;
 485   while (*ptr!='\0' && *ptr<=32) ptr++;
 486   return (*ptr=='\0')? NULL:ptr;
 487 }
 488
 489 // Cuts string at first white space character found by overwriting it with '\0'.
 490 // Returns pointer to next non-white-space char OR to end of string '\0' if none found
 491 inline char* strcut_(char* str)
 492 {
 493   if (!str) return NULL;
 494   char* ptr=str;
 495   while (*ptr!='\0' && *ptr>32) ptr++;
 496   if (*ptr=='\0') return ptr;
 497   *ptr='\0';
 498   ptr++;
 499   while (*ptr!='\0' && *ptr<=32) ptr++;
 500   return ptr;
 501 }
 502
 503 // Cuts string at first occurence of charcter c, by overwriting it with '\0'.
 504 // Returns pointer to next char not equal c, OR to NULL if none found
 505 inline char* strcut(char* str, const char c)
 506 {
 507   if (!str) return NULL;
 508   char* ptr=str;
 509   while (*ptr!='\0' && *ptr!=c) ptr++;
 510   if (*ptr=='\0') return NULL;
 511   *ptr='\0';
 512   ptr++;
 513   while (*ptr!='\0' && *ptr==c) ptr++;
 514   return (*ptr=='\0')? NULL:ptr;
 515 }
 516
 517 // Cuts string at first occurence of charcter c, by overwriting it with '\0'.
 518 // Returns pointer to next char not equal c, OR to end of string '\0' if none found
 519 inline char* strcut_(char* str, const char c)
 520 {
 521   if (!str) return NULL;
 522   char* ptr=str;
 523   while (*ptr!='\0' && *ptr!=c) ptr++;
 524   if (*ptr=='\0') return ptr;
 525   *ptr='\0';
 526   ptr++;
 527   while (*ptr!='\0' && *ptr==c) ptr++;
 528   return ptr;
 529 }
 530
 531 // Cuts string at first occurence of substr, by overwriting the first letter with '\0'.
 532 // Returns pointer to next char after occurence of substr, OR to NULL if no such char found
 533 inline char* strcut(char* str, const char* substr)
 534 {
 535   char* ptr;     //present location in str being compared to substr
 536   const char* sptr=substr; //present location in substr being compared to substr
 537   // while not at end of str and not all of substr is matched yet
 538   while (1)
 539     {
 540       for (ptr=str, sptr=substr; *ptr==*sptr && *ptr!='\0';  ptr++, sptr++) ;
 541       if (*sptr=='\0') {*str='\0'; return ptr;}
 542       if (*ptr=='\0')  return NULL;
 543       str++;
 544     }
 545 }
 546
 547 // Cuts string at first occurence of substr, by overwriting the first letter with '\0'.
 548 // Returns pointer to next char after occurence of substr, OR to end of string '\0' if no such char found
 549 inline char* strcut_(char* str, const char* substr)
 550 {
 551   char* ptr;         //present location in str being compared to substr
 552   const char* sptr=substr; //present location in substr being compared to str
 553   // while not at end of str and not all of substr is matched yet
 554   while (1)
 555     {
 556       for (ptr=str, sptr=substr; *ptr==*sptr && *ptr!='\0';  ptr++, sptr++) ;
 557       if (*sptr=='\0') {*str='\0'; return ptr;}
 558       if (*ptr=='\0')  return ptr;
 559       str++;
 560     }
 561 }
 562
 563 // Copies first word in ptr to str. In other words, copies first block of non whitespace characters,
 564 // beginning at ptr, to str. If a word is found, returns address of second word in ptr or, if no second
 565 // word is found, returns address to end of word ('\0' character) in ptr string. If no word is found
 566 // in ptr NULL is returned.
 567 inline char* strwrd(char* str, char* ptr)
 568 {
 569   ptr=strscn(ptr);    // advance to beginning of next word
 570   if (ptr)
 571     {
 572       while (*ptr!='\0' && *ptr>32) *(str++) = *(ptr++);
 573       *str='\0';
 574       while (*ptr!='\0' && *ptr<=32) ptr++;
 575       return ptr;
 576     }
 577   else return NULL;
 578 }
 579
 580 // Copies first word ***delimited by char c*** in ptr to str. In other words, copies first block of non-c characters,
 581 // beginning at ptr, to str. If a word is found, returns address of second word in ptr or, if no second
 582 // word is found, returns address to end of word ('\0' character) in ptr string. If no word is found
 583 // in ptr NULL is returned.
 584 inline char* strwrd(char* str, char* ptr, const char c)
 585 {
 586   ptr=strscn(ptr,c);    // advance to beginning of next word
 587   if (ptr)
 588     {
 589       while (*ptr!='\0' && *ptr!=c) *(str++) = *(ptr++);
 590       *str='\0';
 591       while (*ptr!='\0' && *ptr==c) ptr++;
 592      return ptr;
 593     }
 594   else return NULL;
 595 }
 596
 597 // Similar to Perl's tr/abc/ABC/: Replaces all chars in str found in one list with characters from the second list
 598 // Returns the number of replaced charactrs
 599 int strtr(char* str, const char oldchars[], const char newchars[])
 600 {
 601   char* ptr;
 602   const char *plist;
 603   int ntr=0;
 604   for (ptr=str; *ptr!='\0'; ptr++)
 605     for (plist=oldchars; *plist!='\0'; plist++)
 606       if (*ptr==*plist)
 607         {
 608           *ptr=newchars[plist-oldchars];
 609           ntr++;
 610           break;
 611         }
 612   return ntr;
 613 }
 614
 615 // Similar to Perl's tr/abc//d: deletes all chars in str found in the list
 616 // Returns number of removed characters
 617 int strtrd(char* str, const char chars[])
 618 {
 619   char* ptr0=str;
 620   char* ptr1=str;
 621   const char *plist;
 622   while (*ptr1!='\0')
 623     {
 624       for (plist=chars; *plist!='\0'; plist++)
 625         if (*ptr1==*plist) break;
 626       if (*plist=='\0') {*ptr0=*ptr1; ptr0++;}
 627       ptr1++;
 628     }
 629   return ptr1-ptr0;
 630 }
 631
 632 // Similar to Perl's tr/a-z//d: deletes all chars in str found in the list
 633 // Returns number of removed characters
 634 int strtrd(char* str, char char1, char char2)
 635 {
 636   char* ptr0=str;
 637   char* ptr1=str;
 638   while (*ptr1!='\0')
 639     {
 640       if (*ptr1>=char1 && *ptr1<=char2) {*ptr0=*ptr1; ptr0++;}
 641       ptr1++;
 642     }
 643   return ptr1-ptr0;
 644 }
 645
 646 // transforms str into an all uppercase string
 647 char* uprstr(char* str)
 648 {
 649   char* s=str;
 650   while (*s !='\0') {if (*s>='a' && *s<='z') *s+='A'-'a';s++;}
 651   return(str);
 652 }
 653
 654 // transforms str into an all uppercase string
 655 char* lwrstr(char* str)
 656 {
 657   char* s=str;
 658   while (*s !='\0') {if (*s>='A' && *s<='Z') *s+='a'-'A'; s++;}
 659   return(str);
 660 }
 661
 662 // transforms chr into an uppercase character
 663 inline char uprchr(char chr)
 664 {
 665   return (chr>='a' && chr<='z')? chr+'A'-'a' : chr;
 666 }
 667
 668 // transforms chr into an lowercase character
 669 inline char lwrchr(char chr)
 670 {
 671   return (chr>='A' && chr<='Z')? chr-'A'+'a' : chr;
 672 }
 673
 674
 675 // Replaces first occurence of str1 by str2 in str. Returns pointer to first occurence or NULL if not found
 676 // ATTENTION: if str2 is longer than str1, allocated memory of str must be long enough!!
 677 inline char* strsubst(char* str, const char str1[], const char str2[])
 678 {
 679   char* ptr = strstr(str,str1);
 680   strcpy(ptr,str2);
 681   return ptr;
 682 }
 683
 684 // Gives elapsed time since first call to this function
 685 inline void ElapsedTimeSinceFirstCall(const char str[])
 686 {
 687   timeval t;
 688   static double tfirst=0;
 689   if (tfirst==0)
 690     {
 691       gettimeofday(&t, NULL);
 692       tfirst = 1E-6*t.tv_usec + t.tv_sec;
 693     }
 694   gettimeofday(&t, NULL);
 695   printf("Elapsed time since first call:%12.3fs %s\n",1E-6*t.tv_usec + t.tv_sec - tfirst,str);
 696 }
 697
 698 // Gives elapsed time since last call to this function
 699 inline void ElapsedTimeSinceLastCall(const char str[])
 700 {
 701   timeval t;
 702   static double tlast=0.0;
 703   if (tlast==0.0)
 704     {
 705       gettimeofday(&t, NULL);
 706       tlast = 1.0E-6*t.tv_usec + t.tv_sec;
 707     }
 708   gettimeofday(&t, NULL);
 709   printf("Elapsed time since last call:%12.3fs %s\n",1.0E-6*t.tv_usec + t.tv_sec - tlast,str);
 710   tlast = 1.0E-6*t.tv_usec + t.tv_sec;
 711 }
 712
 713 inline char* RemovePath(char outname[], char filename[])
 714 {
 715   char* ptr;
 716 #ifdef WINDOWS
 717   ptr=strrchr(filename,92);  //return adress for LAST \ (backslash) in name
 718 #else
 719   ptr=strrchr(filename,'/'); //return adress for LAST / in name
 720 #endif
 721   if (!ptr) ptr=filename; else ptr++;
 722   strcpy(outname,ptr);
 723   return outname;
 724 }
 725
 726 inline char* RemoveExtension(char outname[], char filename[])
 727 {
 728   char *ptr1;
 729   ptr1=strrchr(filename,'.');       //return adress for LAST '.' in name
 730   if (ptr1) {*ptr1='\0'; strcpy(outname,filename); *ptr1='.';} else strcpy(outname,filename);
 731   return outname;
 732 }
 733
 734 inline char* RemovePathAndExtension(char outname[], char filename[])
 735 {
 736   char *ptr, *ptr1;
 737 #ifdef WINDOWS
 738   ptr=strrchr(filename,92);  //return adress for LAST \ (backslash) in name
 739 #else
 740   ptr=strrchr(filename,'/'); //return adress for LAST / in name
 741 #endif
 742   if (!ptr) ptr=filename; else ptr++;
 743   ptr1=strrchr(filename,'.');       //return adress for LAST '.' in name
 744   if (ptr1) {*ptr1='\0'; strcpy(outname,ptr); *ptr1='.';} else strcpy(outname,ptr);
 745   return outname;
 746 }
 747
 748 inline char* Extension(char extension[], char filename[])
 749 {
 750   char* ptr;
 751   ptr=strrchr(filename,'.');      //return adress for LAST '.' in name
 752   if (ptr) strcpy(extension,ptr+1); else *extension='\0';
 753   return extension;
 754 }
 755
 756 // Path includes last '/'
 757 inline char* Pathname(char pathname[], char filename[])
 758 {
 759   char* ptr;
 760   char chr;
 761 #ifdef WINDOWS
 762   ptr=strrchr(filename,92);  //return adress for LAST \ (backslash) in name
 763 #else
 764   ptr=strrchr(filename,'/'); //return adress for LAST / in name
 765 #endif
 766   if (ptr) {chr=*(++ptr); *ptr='\0'; strcpy(pathname,filename); *ptr=chr;} else *pathname='\0';
 767   return pathname;
 768 }
 769
 770 // Swaps two integer elements in array k
 771 inline void swapi(int k[], int i, int j)
 772 {
 773   int temp;
 774   temp=k[i]; k[i]=k[j]; k[j]=temp;
 775 }
 776
 777 // QSort sorting routine. time complexity of O(N ln(N)) on average
 778 // Sorts the index array k between elements i='left' and i='right' in such a way that afterwards
 779 // v[k[i]] is sorted downwards (up=-1) or upwards (up=+1)
 780 void QSortInt(int v[], int k[], int left, int right, int up=+1)
 781 {
 782   int i;
 783   int last;   // last element to have been swapped
 784
 785   if (left>=right) return;        // do nothing if less then 2 elements to sort
 786   // Put pivot element in the middle of the sort range to the side (to position 'left') ...
 787   swapi(k,left,(left+right)/2);
 788   last=left;
 789   // ... and swap all elements i SMALLER than the pivot
 790   // with an element that is LARGER than the pivot (element last+1):
 791   if (up==1)
 792     {
 793       for (i=left+1; i<=right; i++)
 794         if (v[k[i]]<v[k[left]]) swapi(k,++last,i);
 795     }
 796   else
 797     for (i=left+1; i<=right; i++)
 798       if (v[k[i]]>v[k[left]]) swapi(k,++last,i);
 799
 800   // Put the pivot to the right of the elements which are SMALLER, left to elements which are LARGER
 801   swapi(k,left,last);
 802
 803   // Sort the elements left from the pivot and right from the pivot
 804   QSortInt(v,k,left,last-1,up);
 805   QSortInt(v,k,last+1,right,up);
 806 }
 807
 808 // QSort sorting routine. time complexity of O(N ln(N)) on average
 809 // Sorts the index array k between elements i='left' and i='right' in such a way that afterwards
 810 // v[k[i]] is sorted downwards (up=-1) or upwards (up=+1)
 811 void QSortFloat(float v[], int k[], int left, int right, int up=+1)
 812 {
 813   int i;
 814   int last;   // last element to have been swapped
 815   void swapi(int k[], int i, int j);
 816
 817   if (left>=right) return;        // do nothing if less then 2 elements to sort
 818   // Put pivot element in the middle of the sort range to the side (to position 'left') ...
 819   swapi(k,left,(left+right)/2);
 820   last=left;
 821   // ... and swap all elements i SMALLER than the pivot
 822   // with an element that is LARGER than the pivot (element last+1):
 823   if (up==1)
 824     {
 825     for (i=left+1; i<=right; i++)
 826       if (v[k[i]]<v[k[left]]) swapi(k,++last,i);
 827     }
 828   else
 829     for (i=left+1; i<=right; i++)
 830       if (v[k[i]]>v[k[left]]) swapi(k,++last,i);
 831
 832   // Put the pivot to the right of the elements which are SMALLER, left to elements which are LARGER
 833   swapi(k,left,last);
 834
 835   // Sort the elements left from the pivot and right from the pivot
 836   QSortFloat(v,k,left,last-1,up);
 837   QSortFloat(v,k,last+1,right,up);
 838 }
 839
 840 /**
 841  * @brief comparison function for qsort,
 842  * sorts floating point numbers ascendingly
 843  *
 844  * @param cv1 ponter to 1st entry to be sorted
 845  * @param cv2 ponter to 2nd entry to be sorted
 846  *
 847  * @return 0 if entries are equal,
 848  *  +/-1 if 1st greater/smaller than 2nd
 849  */
 850 int CompFltAsc(const void *cv1, const void *cv2){
 851
 852   float f1 = *(float *)cv1;
 853   float f2 = *(float *)cv2;
 854
 855   if      (f1 > f2) { return +1; }
 856   else if (f1 < f2) { return -1; }
 857   else              { return  0; }
 858
 859 } /* this is the end of CompFltAsc() */
 860
 861 //Return random number in the range [0,1]
 862 inline float frand() { return rand()/(RAND_MAX+1.0); }
 863
 864
 865 /////////////////////////////////////////////////////////////////////////////////////
 866 //// Execute system command
 867 /////////////////////////////////////////////////////////////////////////////////////
 868 void runSystem(std::string cmd, int v = 2)
 869 {
 870   if (v>2)
 871     cout << "Command: " << cmd << "!\n";
 872   int res = system(cmd.c_str());
 873   if (res!=0)
 874     {
 875       cerr << endl << "ERROR when executing: " << cmd << "!\n";
 876       exit(1);
 877     }
 878
 879 }