/***************************************************************************
                          VectorManip.cpp  -  description
                             -------------------
    begin                : Don Mai 13 2004
    copyright            : (C) 2004 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 "Types.h"
#include "VectorManip.h"
#include <cmath>
#include <sstream>

template <class T>
void dumpData(const char* text, const vector<T>& vec, ostream& out)
{
   if(!text) return;
   out<<text<<endl;
   for(uint i = 0; i<vec.size(); ++i){
      out<<vec[i]<<" ";
   }
   out<<endl;
}
template <class T>
void dumpDataToBuf(char *buf, size_t cbSize, const char* text, const vector<T>& vec)
{
   if(!buf || !text) return;
   stringstream str(buf, cbSize);
   str<<text<<endl;
   for(uint i = 0; i<vec.size(); ++i){
      str<<vec[i]<<" ";
   }
   str<<endl;

}
template <class T>
void dumpDataToBuf(char *buf, size_t cbSize, const char *text, const char* delim, const vector<T>& vec)
{
   if(!buf || !text) return;
   stringstream str(buf, cbSize);
   str.write(text, strlen(text));
   str.write(delim, strlen(delim));
   for(uint i = 0; i<vec.size(); ++i){
      str<<vec[i]<<delim;
   }
}

template <class T>
void dumpAccents(const char* text, const vector<T>& vec, ostream& out)
{
   if(!text) return;
   out<<text<<endl;
   for(uint i = 0; i<vec.size(); ++i){
      out<<vec[i]<<" ";
      if(i == 7 || i == 11 || i == 15 || i == 18 || i == 22 || i == 24 ) out<<"| ";
   }
   out<<endl;
}


template <class T>
T vectorSum(const vector<T> &vec)
{
   T sum = 0;
   for(uint i = 0; i<vec.size(); ++i){
      sum += vec[i];
   }
   return sum;
}

template <class T>
T vectorWeightedSum(const vector<T> &vec)
{
   T sum = 0;
   for(uint i = 0; i<vec.size(); ++i){
      sum += (vec.size()-i) * vec[i];
   }
   return sum;
}
template <class T>
uint getLast(const vector<T>& vec)
{
   for(uint i = 0; i<vec.size(); ++i){
      if(vec[i] == (T)0) return i;
   }
   return vec.size()-1;
}

template <class T>
vector<T> add(const vector<T> &vec1, const vector<T> &vec2)
{
   vector<T> sum;
   if(vec1.size()!=vec2.size()) return sum;
   for(uint i = 0; i<vec1.size(); ++i){
      sum .push_back(vec1[i]+ vec2[i]);
   }
   return sum;
}

template <class T>
vector<T> sub(const vector<T> &vec1, const vector<T> &vec2)
{
   vector<T> sum;
   if(vec1.size()!=vec2.size()) return sum;
   for(uint i = 0; i<vec1.size(); ++i){
      sum .push_back(vec1[i] - vec2[i]);
   }
   return sum;
}

template <class T>
T scalarProduct(const vector<T> &vec1, const vector<T> &vec2)
{
   T sum = T(0);
   if(vec1.size()!=vec2.size()) return sum;
   for(uint i = 0; i<vec1.size(); ++i){
      sum += vec1[i]*vec2[i];
   }
   return sum;
}

template <class T>
double normEuclid(const vector<T> &vec1)
{
   return sqrt(double(scalarProduct<T>(vec1, vec1)));
}

template <class T>
double cosangle(const vector<T> &vec1, const vector<T> &vec2)
{
   double sc = (double) scalarProduct<T>(vec1, vec2);
   double n1 = normEuclid<T>(vec1);
   double n2 = normEuclid<T>(vec2);
   if(n1 == 0 || n2 == 0) return 0;
   return sc/n1/n2;
}

template <class T>
uint getMaxLen(const vector<T> &vec1, const vector<T> &vec2)
{
   if(vec1.size()<vec2.size()) return vec2.size();
   return vec1.size();
}
template <class T>
uint getMinLen(const vector<T> &vec1, const vector<T> &vec2)
{
   if(vec1.size()>vec2.size()) return vec2.size();
   return vec1.size();
}

template <class T>
vector<T> scaleVector(const vector<T> &vec1, const T &s)
{
   vector<T> sum;
   for(uint i = 0; i<vec1.size(); ++i){
      sum.push_back( vec1[i]* s);
   }
   return sum;
}

template <class T>
vector<T> mult(const vector<T> &vec, const T &s)
{
   vector<T> sum;
   for(uint i = 0; i<vec.size(); ++i){
      sum.push_back( vec[i]* s);
   }
   return sum;
}

template <class T>
vector<T> invScaleVector(const vector<T> &vec1, const T &s)
{
   vector<T> sum;
   for(uint i = 0; i<vec1.size(); ++i){
      sum.push_back( vec1[i] / s);
   }
   return sum;
}

template <class T>
T getMax(const vector<T> &vec, uint& idx)
{
   T max = T(0);
   for(uint i = 0; i<vec.size(); ++i){
      if( i == 0 || vec[i]>max){
         max = vec[i];
         idx = i;
      }
   }
   return max;
}

template <class T>
int find(const vector<T> &vec, const T& val)
{
   for(uint i = 0; i<vec.size(); ++i){
      if(vec[i] == val) return i;
   }
   return -1;
}

template<class T>
vector<T> getSubVector(const vector<T>& vec, uint start, uint end)
{
   vector<T> tmp;
   for(uint i = 0; i<vec.size(); ++i){
      if(i>=start && i<end)
         tmp.push_back(vec[i]);
   }
   return tmp;
}
template<class T>
vector<T> getModus(const vector<T> &vec, uint &modFreq)
{
   uint i = 0, j = 0;
   vector<T> modi;
   vector<pair<T, uint> > freq;
   uint modus = 0;
   modFreq = 1;
   if(vec.empty()) return modi;
   for( i = 0; i<vec.size();++i){
      T val = vec[i];
      //val in modi?
      bool bFound = false;
      for(j= 0; j< freq.size(); ++j){
         if(freq[j].first == val){
            //found: increase frequency
            freq[j].second++;
            if(freq[j].second>modFreq) modFreq = freq[j].second;
            bFound = true;
         }
      }
      //if not found add to list
      if(!bFound){
         pair<T, uint> n;
         n.first = val;
         n.second = 1;
         freq.push_back(n);
      }
   }
   //return vector of values with highest modus
   for(i = 0; i<freq.size(); ++i){
      if(freq[i].second == modFreq)
         modi.push_back(freq[i].first);
   }
   return modi;
}

template<class T>
double getMedian(const vector<T> &vec)
{
   if(vec.empty()) return 0;
   vector<T> tmp = sort<T>(vec, -1);
   double val = 0;
   if(vec.size() % 2){
      val = static_cast<double>(tmp[tmp.size()/2]);
   }
   else{
      val = 0.5*(static_cast<double>(tmp[tmp.size()/2-1]) + static_cast<double>(tmp[tmp.size()/2]));
   }
   return val;
}
template<class T>
uint getNumberOfRuns(const vector<T> vec, const T& med, uint& n1, uint &n2)
{
   unsigned int nRuns = 0; 
   n1 = 0;
   n2 = 0;
   bool bSmallerOld = vec[0]< med ? true : false;
   for(uint i = 1; i<vec.size(); ++i){
      bool bSmaller = vec[i]< med ? true : false;
      if(bSmaller) n2++;
      else n1++;
      if(bSmaller != bSmallerOld){
         nRuns++;
      }
      bSmallerOld = bSmaller;
   }
   nRuns++;
   return nRuns;
}

template<class T>
double getMean(const vector<T>& vec)
{
   double sum = 0;
   uint len = vec.size();
   for(uint i = 0; i< len; ++i){
      sum += static_cast<double>(vec[i]);
   }
   if(len) sum /= (double)len;

   return sum;

}

template<class T>
double getVariance2(const vector<T>& vec)
{
   double sum = 0;
   double mean = getMean<T>(vec);
   uint len = vec.size();
   for(uint i = 0; i< len; ++i){
      sum += (static_cast<double>(vec[i])-mean)*(static_cast<double>(vec[i])-mean);
   }
   if(len) sum = sum / (double)len;
   return sum;

}

template<class T>
double getVariance(const vector<T>& vec)
{
   double sum = 0;
   double mean = 0;
   uint len = vec.size();
   for(uint i = 0; i< len; ++i){
      sum += static_cast<double>(vec[i])*static_cast<double>(vec[i]);
      mean += static_cast<double>(vec[i]);
   }
   if(len){
       sum  /= (double)len;
       mean /= (double)len;
   }
   return sum-mean*mean;

}

template<class T>
double getSkewness(const vector<T> &vec)
{
   double sum = 0;
   uint len = vec.size();
   if(len == 0) return 0;
   double mean = getMean<T>(vec);
   double var = getVariance<T>(vec);
   for(uint i = 0; i< len; ++i){
      double v  = static_cast<double>(vec[i]);
      sum += (v-mean)*(v-mean)*(v-mean);
   }

   var = pow(var, 1.5);
   double num = sum/len;
   return num/var;
}

template<class T>
double getKurtosis(const vector<T> &vec)
{
   double sum = 0;
   uint len = vec.size();
   if(len == 0) return 0;
   double mean = getMean<T>(vec);
   double var = getVariance<T>(vec);
   for(uint i = 0; i< len; ++i){
      double v  = static_cast<double>(vec[i]);
      sum += (v-mean)*(v-mean)*(v-mean)*(v-mean);
   }
   var = var *var;
   double num = sum/len;
   return num/var;
}

template<class T>
vector< pair<T, int> > getDistribution(const vector<T> &vec)
{
   vector<std::pair<T, int> > ret;
   
   double sum = 0;
   bool bFound = false;
   
   uint len = vec.size();
   if(len == 0) return ret;

   for(uint i = 0; i< len; ++i){
      bool bFound = false;
      //cerr<<"Val: "<<vec[i];

      for(uint j = 0; j< ret.size(); ++j){
         if(ret[j].first == vec[i]) {
            ret[j].second++;
            bFound = true;
            //cerr<<" found. New freq "<<ret[j].second<<endl;
            break;
         } 
      }
      if(!bFound){
         //cerr<<" added."<<endl;
         ret.push_back(pair<T, int>(vec[i], 1));
      }
   }   
   return ret;
}

template<class T>
double getEntropy(const vector<T> &vec, bool bNorm = true)
{
   double sum = 0;
   uint len = vec.size();
   if(len == 0) return 0;
   vector<pair<T, int> > dist = getDistribution(vec);
   for(uint i = 0; i< dist.size(); ++i){
      //cerr<<"Val: "<<dist[i].first<<" Freq="<<dist[i].second<<endl;
      double p = (dist[i].second +.0)/ len;
      double entropy = p * log((double)p) / log(2.0);
      sum += entropy;
   }
   //double norm = log(dist.size()+.0)/log(2.0);
   double norm = bNorm? log(dist.size()+.0)/log(2.0) : 1;
   return norm? -sum/norm : -sum;
}



template<class T>
void scale(vector<T>& vec, const T scale)
{
   uint i = 0;
   for(i = 0; i<vec.size(); ++i){
      vec[i] *= scale;
   }
}
template<class T>
void invScale(vector<T>& vec, const T scale)
{
   uint i = 0;
   for(i = 0; i<vec.size(); ++i){
      vec[i] /= scale;
   }
}


template<class T>
T getMaximum(const vector<T>& vec)
{
   uint dummy;
   return getMaximum<T>(vec, dummy);
}

template<class T>
T getMaximum(const vector<T>& vec, uint& pos)
{
   uint i = 0;

   if(vec.empty()) return T(0);

   T max = vec[0];

   for(i = 0; i<vec.size();++i){
      if(vec[i]>max){
          max = vec[i];
          pos = i;
      }
   }

   return max;
}

template<class T>
T getMinimum(const vector<T>& vec)
{
   uint dummy;
   return getMinimum<T>(vec, dummy);
}

template<class T>
T getMinimum(const vector<T>& vec, uint& pos)
{
   uint i = 0;
   if(vec.empty()) return T(0);
   T min = vec[0];

   for(i = 0; i < vec.size();++i){
      if(vec[i]<min){
          min = vec[i];
          pos = i;
      }
   }
   return min;
}
template<class T>
vector<T>  getRow(vector<vector<T> > vecvec, uint row)
{
   vector<T> ret;
   if(row>vecvec.size()) return ret;
   return vecvec[row];
   
}

template<class T>
vector<T>  getColumn(vector<vector<T> > vecvec, uint column)
{
   vector<T> ret;
   uint i = 0;
   for(uint i = 0; i<vecvec.size(); ++i){
      if(column<vecvec[i].size()){
         ret.push_back((vecvec[i])[column]);
      }
      else{
         ret.push_back((T)0);       
      }
   }
   return ret;
}


template<class T>
vector<T> getClasses(const vector<T>& vec, const vector<T>& classBorders)
{
   uint i = 0, j = 0, max = classBorders.size()-1;
   vector<T> ret;
   
   for(i = 0; i<vec.size();++i){
      if(vec[i]>=classBorders[max]) {
         ret.push_back(max);
         continue;
      }
      for(j = 0; j<max; ++j){
         if(vec[i]>=classBorders[j] && vec[i]<classBorders[j+1]){
            ret.push_back(j);
            break;
         }
         
      }
   }
   return ret;
}

template<typename T>
vector<T>& remove(vector<T>& vec, const uint idx)
{
   if(idx>=vec.size()) return vec;
   typename vector<T>::iterator it;
   uint i = 0;
   for(it  = vec.begin() ; it != vec.end(); ++it, ++i){
      if(i == idx) {
         vec.erase(it);
         break;
      }
   }
   return vec;
}

static uint getRandomNumber(const uint max)
{
   return (uint)floor( max * ( (rand() +0.0)/RAND_MAX ));
}


template<class T>
vector<T> getRandomSelection(const vector<T>& vec, const uint count)
{
   vector<T> ret;
   vector<T> dummy = vec;

   if(count>vec.size()) return vec;

   for(uint i = 0; i<count; ++i){
      uint idx = getRandomNumber(dummy.size()-1);
      //cout<<"RNo = "<<idx<<endl;
      ret.push_back(dummy[idx]);
      remove(dummy, idx);      
   }
   cout<<endl;
   return ret;
}
template<class T>
bool isEqual(const vector<T>& data1,const vector<T>& data2)
{
   if(data1.size() != data2.size()) return false;
   for(uint i = 0; i<data1.size(); ++i){
      if(data1[i] != data2[i]) return false;
   }
   return true;
}

template<class T>
uint getElementCount(const vector<T>& data, const T& lowerBound, const T& upperBound, int dir = 0)
{
   uint nCount = 0;
   for(size_t i = 0; i< data.size(); ++i){
      if(dir >0 && data[i]<upperBound && data[i]>lowerBound
         || 
         dir <0 && (data[i]>upperBound || data[i]<lowerBound)         
         ) nCount++;
   }
   return nCount;
}

template<class T>
uint getElementCount(const vector<T>& data, const T& threshold, int dir)
{
   uint nCount = 0;
   for(size_t i = 0; i< data.size(); ++i){
      if(dir > 0 && data[i]>=threshold 
         ||
         dir < 0 && data[i]<=threshold 
         ||
         dir == 0 && data[i]==threshold 
         ) nCount++;
   }
   return nCount;

}

template<class T>
bool Greater(const T& e1,const T& e2)
{
   if(e1<e2) return false;
   return true;
} 


template<class T>
vector<T> sort(const vector<T>& data, int dir)
{
   vector<T> ret(data);
   uint s =ret.size();
   if(dir >0) sort(ret.begin(), ret.end(), Greater<T>);
   else std::sort(ret.begin(), ret.end());
   //dumpData("Sorted:", ret,cout);
   return ret;
}
template<class T>
vector<T> fuzzyNot(const vector<T>& vec)
{
   vector<T> ret;
   for(uint i = 0; i<vec.size(); ++i){
      ret.push_back(T(1) - vec[i]);
   }
   return ret;
}

template<class T>
vector<T> fuzzyOr(const vector<T>& vec1, const vector<T>& vec2)
{
   vector<T> ret;
   for(uint i = 0; i< vec1.size(); ++i){
      T val;
      if(i>vec2.size()-1) val = 0;
      else val = vec2[i];
      ret.push_back(max<T>(vec1[i], val));
   }
   return ret;
}

template<class T>
vector<T> fuzzyAnd(const vector<T>& vec1, const vector<T>& vec2)
{
   vector<T> ret;
   for(uint i = 0; i< vec1.size(); ++i){
      T val;
      if(i>vec2.size()-1) val = 0;
      else val = vec2[i];
      ret.push_back(min<T>(vec1[i], val));
   }
   return ret;
}

template<class T>
double  match(const vector<T>& vec1, const vector<T>& vec2, bool bNorm)
{
   double ret = -1.0;
   if(vec1.size() != vec2.size() || vec1.empty()) return ret;
   uint matches = 0;
   for(uint i = 0 ; i<vec1.size(); ++i){
      if(vec1[i] == vec2[i]) matches++;
   }
   return bNorm ? double(matches)/double(vec1.size()) : double(matches);
}

template<class T>
double  mismatch(const vector<T>& vec1, const vector<T>& vec2, bool bNorm )
{
   double ret = -1.0;
   if(vec1.size() != vec2.size() || vec1.empty()) return ret;
   uint mismatches = 0;
   for(uint i = 0 ; i<vec1.size(); ++i){
      if(vec1[i] != vec2[i]) mismatches++;
   }
   return bNorm ? double(mismatches)/double(vec1.size()) : double(mismatches);

}
