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