/***************************************************************************
                          Misc.cpp  -  description
                             -------------------
    begin                : Sam Feb 8 2003
    copyright            : (C) 2003 by klaus frieler
    email                : kgf@omniversum.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#include <cmath>
#include "Misc.h"
#include "gauss.h"

double  getSum(const data& vec)
{
   const_dataIter it;
   double sum = 0.0;

   for(it = vec.begin(); it != vec.end();++it){
      sum += *it;
   }

   return sum;
}

double  getSquareSum(const data& vec)
{
   const_dataIter it;
   double sum = 0.0;

   for(it = vec.begin(); it != vec.end();++it){
      sum += (*it) * (*it);
   }

   return sum;
}

double getMinimum(const data& vec)
{
   uint dummy;
   return getMinimum(vec, dummy);
}

double getMinimum(const data& vec, uint& pos)
{
   const_dataIter it;
   double min = *vec.begin();

   for(it = vec.begin(); it != vec.end();++it){
      if(*it<min){
          min = *it;
          pos = (uint) ( it - vec.begin());
      }
   }
   return min;
}

double getMaximum(const data& vec)
{
   uint dummy;
   return getMaximum(vec, dummy);
}

double getMaximum(const data& vec, uint& pos)
{
   
   const_dataIter it;
   double max = *vec.begin();

   for(it = vec.begin(); it != vec.end();++it){
      if(*it>max){
          max = *it;
          pos = (uint) (it-vec.begin());
      }
   }
   
   return max;
}

uint getMinimumN(double est, double mean, double thres, uint max)
{
   uint i = max + 1 ;
   for ( i = 1; i<max; ++i){
      if( fabs(i*est - mean)<thres) break;
   }

   return i;
}

double calcSquareErr(double est, double thres, const data& vec, uint max, double& N )
{
   const_dataIter it;
   data           idx;
   
   double   sum  = 0;
   double   h = 0;
   uint     i = 0;
   
   for(it = vec.begin();it!=vec.end(); ++it){
      i = getMinimumN(est, *it, thres, max);
      if ( i == max + 1) return -1.0;
      idx.push_back((double)i);
   }

   for(i = 0; i<idx.size();++i){
     h    = idx[i]*est-vec[i];
     sum += h*h;
   }

   N = getSum(idx);

   return sum;
}


void dumpData(const data& d, const char *delim)
{
   for(uint i = 0; i<d.size(); ++i){
      cout<<d[i]<<delim;
   }
   cout<<endl;
}

void dumpData(const vector<uint>& d, const char* delim)
{
   for(uint i = 0; i<d.size(); ++i){
      cout<<d[i]<<delim;
   }
   cout<<endl;
}
void dumpData(const data& d, const vector<uint>& p, const char* delim)
{
   if(d.size() != p.size()) return;
   for(uint i = 0; i<d.size(); ++i){
   cout<<p[i]<<delim<<d[i]<<endl;
   }
}

double calcMinimumDistance(const data &vec, uint &i, uint &j)
{
   double diff = 0, minDiff = 0;
   uint minI = 0,  minJ = 0;
   i = 0;
   j = 0;
   if (vec.empty()) return false;
   minDiff = fabs(vec[0]-vec[1]);
   for (uint k = 0; k<vec.size(); ++k)
      for (uint l = k+1; l<vec.size(); ++l){
         diff = fabs(vec[k]-vec[l]);
         if (diff<=minDiff){
            minDiff = diff;
            minI = k;
            minJ = l;
         }
      }

   i = minI;
   j = minJ;
      
   return minDiff;
}


char *diffString(const char *str1, const char *str2)
{
   static char buf[1024];
   uint len = strlen(str1);
   if ( strlen(str2)<len) len = strlen(str2);  

   if (len>sizeof(buf)-1) return 0;

   for (uint i = 0; i<len;++i){
      if( str1[i] == str2[i])
         buf[i] = '0';
      else
         buf[i] = '1';
   }

   buf[len] = 0;
   return buf;
}

double match(const double *vec1, const double *vec2, const uint len)
{
   double sum = 0.0;
   double n1  = norm(vec1,len), n2 = norm(vec2,len);


   for (uint i = 0;i<len;i++){
      sum += vec1[i] * vec2[i];
   }

   if (n1 && n2)
      return sum/sqrt(n1*n2);

   return 0.0;
}

double match2(const double *vec1, const double *vec2, const uint len)
{
   double sum = 0.0;
   double n1  = 0.0, n2 = 0.0;

   for (uint i = 0; i<len; i++){
      sum += vec1[i] * vec2[i] * exp(-(i+0.0)/len);
      n1  += vec1[i] * vec1[i] * exp(-(i+0.0)/len);
      n2  += vec2[i] * vec2[i] * exp(-(i+0.0)/len);
   }


   if (n1 && n2)
      return sum/sqrt(n1*n2);

   return 0.0;
}

double match2(const double *vec1, const double val, const uint len)
{
   double sum = 0.0;
   double n1  = 0.0, n2 = 0.0;

   for (uint i = 0; i<len; i++){
      sum += vec1[i] * val     * exp(-(i+0.0)/len);
      n1  += vec1[i] * vec1[i] * exp(-(i+0.0)/len);
      n2  += val     * val     * exp(-(i+0.0)/len);
    }


   if (n1 && n2)
      return sum/sqrt(n1*n2);

   return 0.0;
}

double mismatch2(const double *vec, const double val, const uint len)
{
   double sum = 0.0;
   for (uint i = 0; i<len; i++)
      sum += (vec[i] - val) * (vec[i] - val);

   sum = sqrt(sum/(val * val *(double) len));

   return sum;
}


double mismatch3(const double *vec1, const double *vec2, const uint len)
{
   double sum = 0.0;
   double n1 = norm(vec1, len), n2 = norm(vec2, len);
   
   for (uint i = 0; i<len; i++){
      double tmp = (vec1[i]/n1 - vec2[i]/n2);
      sum += tmp * tmp;
   }

   sum = sqrt(sum);

   return sum;
}

double norm(const double *vec,  const uint len)
{
   double sum = 0.0;
   for (uint i = 0; i<len; i++)
      sum += (vec[i] * vec[i]);

   return sqrt(sum);
}

void normalize(double *vec, const uint len)
{
   double n = norm(vec, len);

   if (n == 0) return;
   for (uint i = 0; i<len; i++)
      vec[i] = vec[i]/n;

   return;
}


bool  writeDataToFile(string file, const data& vec, double step, double start)
{
   bool bRet = false;
   FILE *fp;
   const_dataIter it;
   double  t = start;
   if( NULL == ( fp = fopen(file.c_str(),"w")))
      return bRet;

   it = vec.begin();
   while( !feof(fp) && it != vec.end()){
      char buf[128];
      sprintf(buf,"%f %f\n", t, *it++);
      fputs(buf,fp);
      if(ferror(fp)) goto out;
      t += step;
   }
   bRet = true;

out:
   fclose(fp);
   return bRet;  
}

bool writeDataToFile(string file, const data& vec,  const data& step)
{
   bool bRet = false;
   FILE *fp;
   const_dataIter it;
   const_dataIter it1;

   if( NULL == ( fp = fopen(file.c_str(),"w")))
      return bRet;

   it    = vec.begin();
   it1   = step.begin();
   
   while( !feof(fp) && it != vec.end() && it1 != step.end()){
      char buf[128];
      sprintf(buf,"%f %f\n", *it1++, *it++);
      fputs(buf,fp);
      if(ferror(fp)) goto out;
   }
   bRet = true;

out:
   fclose(fp);
   return bRet;
}

void scale(data& vec, double scale)
{
   dataIter it;
   for(it = vec.begin(); it != vec.end(); ++it){
      *it *= scale;
   }
}

data scale(const vector<uint>& vec, double scale)
{
   data tmp;
   vector<uint>::const_iterator it;
   for(it = vec.begin(); it != vec.end(); ++it){
      tmp.push_back(*it * scale);
   }
   return tmp;
}

data getMaxima(const data& vec, vector<uint>& points, data& values)
{
   data tmp;
   uint size = vec.size();
   values.clear();
   points.clear();
   if(vec[0]>=vec[1]){
      points.push_back(0);
      values.push_back(vec[0]);
   }
   for (uint i = 1; i <size-2;++i){
      if(vec[i]>vec[i-1] && vec[i]>vec[i+1]){
         points.push_back(i);
         values.push_back(vec[i]);
      }
   }

   if(vec[size-1]>vec[size-2]){
      points.push_back(size-1);
      values.push_back(vec[size-1]);
   }
   return values;
}

data getMinima(const data& vec, vector<uint>& points, data& values)
{
   data tmp;
   uint size = vec.size();
   values.clear();
   points.clear();

   if(vec[0]<vec[1]){
      points.push_back(0);
      values.push_back(vec[0]);
   }
   for (uint i = 1; i<size-2; ++i){
      if(vec[i]<vec[i-1] && vec[i]<vec[i+1]){
         points.push_back(i);
         values.push_back(vec[i]);
      }
   }

   if(vec[size-1]>vec[size-2]){
      points.push_back(size-1);
      values.push_back(vec[size-1]);
   }

   return values;
}

void getMaxPair(const data& values, const vector<uint>& points, double &val, uint &pos )
{
   uint pos1 = 0;
   val = getMaximum(values, pos1);
   pos = points[pos];
}


void   calcIOIs(const data& timePoints, data& iois)
{
   const_dataIter it;
   iois.clear();

   for(it = timePoints.begin()+1; it != timePoints.end(); ++it){
      iois.push_back(*it-*(it-1));
   }
}

void calcTimePoints(const data& iois, data& timePoints)
{
   const_dataIter it;
   timePoints.clear();
   double time = 0;
   timePoints.push_back(time);
   for(it = iois.begin(); it != iois.end(); ++it){
      time += *it;
      timePoints.push_back(time);
   }
}

void shiftData(data & vec, double offset)
{
   
   for(dataIter it = vec.begin(); it != vec.end(); ++it){
      *it -= offset;
   }

}

bool isNumber(const char *str)
{
   const char *numeric = "0123456789.";
   if(!str) return false;
   if(str[0] == '\n' || str[0] == '\r') return false;
   for (uint i = 0; i<strlen(str)-1; ++i){
      if(!strchr(numeric, (int) str[i]))return false;      
   }
   return true;
}


void  getAutoCorrMaxima(const data &autoCorr, data& dMaxPoints, data& maxima, double autoCorrStep, const bool verbose, const bool correct)
{
   data maxima_corr;
   vector<uint> maxPoints;

   getMaxima(autoCorr, maxPoints, maxima);
   //dumpData(maxPoints);
   //dumpData(maxima);
   dMaxPoints = scale( maxPoints, autoCorrStep);
   //dumpData(dMaxPoints);
   if(correct){
      double dx = dMaxPoints[1]-dMaxPoints[0];
      double dy = -1.0/(autoCorr.size());
      for (uint i = 0; i<maxima.size(); ++i){
         double e = dy/dx*dMaxPoints[i] + 1;
         maxima_corr.push_back(maxima[i] / e);
      }
      maxima = maxima_corr;
   }

   if(verbose){
      cout<<"Gauss auto correlation:"<<endl;
      cout<<"Point   Value   Corrected value"<<endl;
      cout<<"-------------------------------"<<endl;
      for(uint i = 0; i<maxima.size(); ++i){
         char buf[256];
         if(correct){
            sprintf(buf, "%05.f   %1.2f    %1.2f",dMaxPoints[i], maxima[i], maxima_corr[i]);
         }
         else{
            sprintf(buf, "%05.f   %1.2f    ---",dMaxPoints[i], maxima[i]);
         }
         cout<<buf<<endl;
      }
      writeDataToFile(string("auto_max.dat"), maxima , dMaxPoints);

   }
}

void evalAutoCorrMaxima(const data &max, const data& dMaxPoints, uint &tactus, uint &period)
{
   data maxima(max);
   maxima[0]  = 0;//alway max in auto corr, clear it
   uint pos = 0;
   getMaximum(maxima, pos);
   tactus = (int)floor(0.5+dMaxPoints[pos]/dMaxPoints[1]);

   period = 0;
   while(period<tactus){
      maxima[pos] = 0;//for getting second max
      getMaximum(maxima, pos);
      period = (int)floor(0.5+dMaxPoints[pos]/dMaxPoints[1]);
   }

}

string toLower(const string& str)
{
   string sRet("");
   for(uint i = 0; i<str.length(); ++i){
      if( str[i] >= 'A' && str[i]<='Z'){
         sRet.push_back(str[i]-'A' +'a');
      }
      else sRet.push_back(str[i]);
   }
   return sRet;
}

string   stripWhiteSpace(const string& str, bool left, bool right)
{
   uint idx = 0;
   string tmp(str);
   if(left){
      idx = tmp.find_first_not_of(" \n\t\r");
      if( idx != string::npos){
         tmp = tmp.substr(idx);
      }
   }
   if(right){
      idx = tmp.find_last_not_of(" \t\n\r");
      if( idx != string::npos){
         tmp = tmp.substr(0,idx + 1);
      }
      
   }
   return tmp;
}

string normalizeString(const string &str)
{
   return stripWhiteSpace(toLower(str), true, true);
}

string getToken(const char *str, uint which, const char delim)
{
   uint i = 0;
   char *tmp;
   char *start = 0, *end = 0;
   string szRet("");
   if(!str || !which) return szRet;

   tmp = new char[strlen(str)+1];
   if(!tmp) return szRet;
   strcpy(tmp, str);

   start = tmp;
   end   = strchr(start, delim);

   if(!end){
      if(which != 1){
         goto out;
      }
      else{
         szRet = str;
         goto out;
      }
   }

   for(i = 0; i<which - 1; ++i){
      end = strchr(start, delim);
      if(!end){
         if(i == which-1){
            end = start + strlen(start);
            break;
         }
         else goto out;
      }
      start = end+1;
   }
   end = strchr(start, delim);
   if(end)*end = '\0';
   szRet = string(start);

out:
   delete []tmp;
   return szRet;
}

string stripPath(const char* fullPath)
{
   if(!fullPath) return string("");
   string tmp = fullPath;
   int pos = tmp.find_last_of('\\');
   if(pos == -1) return tmp;
   return tmp.substr(pos+1);
}
string stripExtension(const char* name)
{
   if(!name) return string("");
   string tmp = name;
   int pos = tmp.find_last_of('.');
   if(pos == -1) return tmp;
   return tmp.substr(0, pos);   
}
char *subchr(char *str, char from, char to)
{
   if(!str) return str;
   uint len = strlen(str);
   
   for(uint i = 0; i<len; ++i){
      if(str[i] == from)
         str[i] = to;
   }
   return str;
}

string   substituteChar(const string& str, const char which, const char what)
{
   string buf(str);
   for(size_t i = 0; i<buf.length(); ++i){
      if(buf[i] == which) buf[i] = what;
   }
   return buf;
}

bool SimType::operator ==(const SimType& comp)
{ 
  return toLower(m_sName) == toLower(comp.m_sName);
}
bool  SimType::operator!=(const SimType& comp){
  return ! (this->operator==(comp));
}

uint ggT(uint m, uint n)
{
  uint z = 0;
  while (n>0){
    z = n;
    n = m % n;
    m = z;
  }
  return m;
}

uint kgV(uint m, uint n)
{
  return (m*n)/ggT(m,n);
}

vector<uint> getPrimeFactorDecomposition(const uint n)
{
   uint max = (uint)floor(sqrt((double)n)) + 1;
   uint p = 2;
   vector<uint> ret;
   
   if(n == 2 || n == 3){
      ret.push_back(n);
      return ret;
   }
   
   while(n % p ){
      p++;
      if(p >= max){ 
         ret.push_back(n);
         return ret;
      }
   }
   vector<uint> g = getPrimeFactorDecomposition(p);
   ret.insert(ret.end(), g.begin(), g.end());
   g = getPrimeFactorDecomposition(n/p);
   ret.insert(ret.end(), g.begin(), g.end());
   return ret;
}
vector<char>   stringToVector(const string &str)
{
   vector<char> ret;
   for(uint i = 0; i<str.size(); ++i){
      ret.push_back(str[i]);
   }
   return ret;
}
