VectorTools.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1

#ifndef BPP_NUMERIC_VECTORTOOLS_H
#define BPP_NUMERIC_VECTORTOOLS_H


#include "../App/ApplicationTools.h"
#include "../Io/OutputStream.h"
#include "AdaptiveKernelDensityEstimation.h"
#include "Matrix/Matrix.h"
#include "NumTools.h"
#include "VectorExceptions.h"

// From the STL:
#include <vector>
#include <map>
#include <cmath>
#include <algorithm>
#include <complex>
#include <numeric>

namespace bpp
{
typedef std::vector<std::complex<double>> Vcomplex;
typedef std::vector<Vcomplex> VVcomplex;
typedef std::vector<VVcomplex> VVVcomplex;

typedef std::vector<std::complex<long double>> Vlcomplex;
typedef std::vector<Vlcomplex> VVlcomplex;
typedef std::vector<VVlcomplex> VVVlcomplex;

typedef std::vector<double> Vdouble;
typedef std::vector<Vdouble> VVdouble;
typedef std::vector<VVdouble> VVVdouble;
typedef std::vector<VVVdouble> VVVVdouble;

typedef std::vector<long double> Vldouble;
typedef std::vector<Vldouble> VVldouble;
typedef std::vector<VVldouble> VVVldouble;
typedef std::vector<VVVldouble> VVVVldouble;

typedef std::vector<int> Vint;
typedef std::vector<Vint> VVint;
typedef std::vector<VVint> VVVint;
typedef std::vector<VVVint> VVVVint;

typedef std::vector<unsigned int> Vuint;
typedef std::vector<Vuint> VVuint;
typedef std::vector<VVuint> VVVuint;
typedef std::vector<VVVuint> VVVVuint;

template<class T>
std::vector<T>  operator+(const std::vector<T>& v1, const std::vector<T>& v2)
{
  size_t size;
  if (v1.size() != v2.size())
  {
    throw DimensionException("VectorTools::operator+", v1.size(), v2.size());
  }
  else
  {
    size = v1.size();
  }
  std::vector<T> result(size);
  for (size_t i = 0; i < size; i++)
  {
    result[i] = v1[i] + v2[i];
  }
  return result;
}

template<class T>
std::vector<T> operator-(const std::vector<T>& v1, const std::vector<T>& v2)
{
  size_t size;
  if (v1.size() != v2.size())
  {
    throw DimensionException("VectorTools::operator-", v1.size(), v2.size());
  }
  else
  {
    size = v1.size();
  }
  std::vector<T> result(size);
  for (size_t i = 0; i < size; i++)
  {
    result[i] = v1[i] - v2[i];
  }
  return result;
}

template<class T>
std::vector<T> operator*(const std::vector<T>& v1, const std::vector<T>& v2)
{
  size_t size;
  if (v1.size() != v2.size())
  {
    throw DimensionException("VectorTools::operator*", v1.size(), v2.size());
  }
  else
  {
    size = v1.size();
  }
  std::vector<T> result(size);
  for (size_t i = 0; i < size; i++)
  {
    result[i] = v1[i] * v2[i];
  }
  return result;
}

template<class T>
std::vector<T> operator/(const std::vector<T>& v1, const std::vector<T>& v2)
{
  size_t size;
  if (v1.size() != v2.size())
  {
    throw DimensionException("VectorTools::operator/", v1.size(), v2.size());
  }
  else
  {
    size = v1.size();
  }
  std::vector<T> result(size);
  for (size_t i = 0; i < size; i++)
  {
    result[i] = v1[i] / v2[i];
  }
  return result;
}


template<class T, class C>
std::vector<T> operator+(const std::vector<T>& v1, const C& c)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = v1[i] + T(c);
  }
  return result;
}
template<class T, class C>
std::vector<T> operator+(const C& c, const std::vector<T>& v1)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = T(c) + v1[i];
  }
  return result;
}

template<class T, class C>
std::vector<T> operator-(const std::vector<T>& v1, const C& c)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = v1[i] - T(c);
  }
  return result;
}
template<class T, class C>
std::vector<T> operator-(const C& c, const std::vector<T>& v1)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = T(c) - v1[i];
  }
  return result;
}

template<class T, class C>
std::vector<T> operator*(const std::vector<T>& v1, const C& c)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = v1[i] * c;
  }
  return result;
}
template<class T, class C>
std::vector<T> operator*(const C& c, const std::vector<T>& v1)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = c * v1[i];
  }
  return result;
}

template<class T, class C>
std::vector<T> operator/(const std::vector<T>& v1, const C& c)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = v1[i] / c;
  }
  return result;
}
template<class T, class C>
std::vector<T> operator/(const C& c, const std::vector<T>& v1)
{
  std::vector<T> result(v1.size());
  for (size_t i = 0; i < result.size(); i++)
  {
    result[i] = c / v1[i];
  }
  return result;
}

template<class T>
void operator+=(std::vector<T>& v1, const std::vector<T>& v2)
{
  for (size_t i = 0; i < v1.size(); i++)
  {
    v1[i] += v2[i];
  }
}

template<class T>
void operator-=(std::vector<T>& v1, const std::vector<T>& v2)
{
  for (size_t i = 0; i < v1.size(); i++)
  {
    v1[i] -= v2[i];
  }
}

template<class T>
void operator*=(std::vector<T>& v1, const std::vector<T>& v2)
{
  for (size_t i = 0; i < v1.size(); i++)
  {
    v1[i] *= v2[i];
  }
}

template<class T>
void operator/=(std::vector<T>& v1, const std::vector<T>& v2)
{
  for (size_t i = 0; i < v1.size(); i++)
  {
    v1[i] /= v2[i];
  }
}

template<class T, class C>
void operator&=(std::vector<T>& v1, const C& c)
{
  for (auto& x : v1)
  {
    x = c;
  }
}

template<class T, class C>
void operator+=(std::vector<T>& v1, const C& c)
{
  for (auto& x : v1)
  {
    x += c;
  }
}

template<class T, class C>
void operator-=(std::vector<T>& v1, const C& c)
{
  for (auto& x : v1)
  {
    x -= c;
  }
}

template<class T, class C>
void operator*=(std::vector<T>& v1, const C& c)
{
  for (auto& x :v1)
  {
    x *= c;
  }
}

template<class T, class C>
void operator/=(std::vector<T>& v1, const C& c)
{
  for (auto& x : v1)
  {
    x /= c;
  }
}

/******************************************************************************/

class VectorTools
{
public:
  VectorTools() {}
  virtual ~VectorTools() {}

public:
  static void resize2(VVdouble& vv, size_t n1, size_t n2)
  {
    vv.resize(n1);
    for (auto& x : vv) { x.resize(n2); }
  }

  static void resize3(VVVdouble& vvv, size_t n1, size_t n2, size_t n3)
  {
    vvv.resize(n1);
    for (auto& vv : vvv)
    {
      vv.resize(n2);
      for (auto& v : vv)
      {
        v.resize(n3);
      }
    }
  }

  static void resize4(VVVVdouble& vvvv, size_t n1, size_t n2, size_t n3, size_t n4)
  {
    vvvv.resize(n1);
    for (auto& vvv : vvvv)
    {
      vvv.resize(n2);
      for (auto& vv : vvv)
      {
        vv.resize(n3);
        for (auto& v : vv)
        {
          v.resize(n4);
        }
      }
    }
  }

  template<class T>
  static void fill(std::vector<T>& v, T value)
  {
    std::fill(v.begin(), v.end(), value);
  }

  template<class T>
  static std::vector<T> seq(T from, T to, T by)
  {
    std::vector<T> v((size_t)((std::abs(from - to) + by / 100) / by) + 1);
    T step = from < to ? by : -by;
    T val(from < to ? from : to);
    for (auto& vi:v)
    {
      vi = val;
      val += step;
    }

    return v;
  }

  template<class T>
  static size_t which(const std::vector<T>& v, const T& which)
  {
    for (size_t i = 0; i < v.size(); i++)
    {
      if (v[i] == which) return i;
    }
    throw ElementNotFoundException<T>("VectorTools::which.", &v, &which);
  }

  template<class T>
  static std::vector<size_t> whichAll(const std::vector<T>& v, const T& which)
  {
    std::vector<size_t> w;
    for (size_t i = 0; i < v.size(); i++)
    {
      if (v[i] == which) w.push_back(i);
    }
    if (w.size())
      return w;
    throw ElementNotFoundException<T>("VectorTools::whichAll.", &v, &which);
  }

  template<class T>
  static std::vector<T> unique(const std::vector<T>& v)
  {
    if (v.size() == 0) return v;
    std::vector<T> sortedV(v.begin(), v.end());
    sort(sortedV.begin(), sortedV.end());
    std::vector<T> uniq;
    uniq.push_back(sortedV[0]);
    for (size_t i = 1; i < sortedV.size(); i++)
    {
      if (sortedV[i] != sortedV[i - 1]) uniq.push_back(sortedV[i]);
    }
    return uniq;
  }

  template<class T>
  static bool isUnique(const std::vector<T>& v)
  {
    if (v.size() == 0) return true;
    std::vector<T> sortedV(v.begin(), v.end());
    sort(sortedV.begin(), sortedV.end());
    for (size_t i = 1; i < sortedV.size(); i++)
    {
      if (sortedV[i] == sortedV[i - 1]) return false;
    }
    return true;
  }

  template<class T>
  static std::vector<T> extract(const std::vector<T>& v1, const std::vector<size_t>& v2)
  {
    std::vector<T> v(v2.size());
    for (size_t i = 0; i < v2.size(); ++i)
    {
      v[i] = v1[v2[i]];
    }
    return v;
  }

  template<class T>
  static std::map<T, size_t> countValues(const std::vector<T>& v)
  {
    std::map<T, size_t> c;
    for (auto x : v)
    {
      c[x]++;
    }
    return c;
  }

  static std::vector<double> breaks(const std::vector<double>& v, unsigned int n);

  template<class T>
  static size_t nclassScott(const std::vector<T>& v)
  {
    std::vector<T> r1 = VectorTools::range(v);
    T r = r1[1] - r1[0];
    double n = v.size();
    double h = 3.5 * VectorTools::sd<T, double>(v) * std::pow(n, -1. / 3);
    return (size_t) ceil(r / h);
  }

  template<class T>
  static T prod(const std::vector<T>& v1)
  {
    T p = 1;
    for (auto x : v1)
    {
      p *= x;
    }
    return p;
  }

  template<class T>
  static std::vector<T> cumProd(const std::vector<T>& v1)
  {
    std::vector<T> p(v1.size());
    if (v1.size() == 0) return p;
    p[0] = v1[0];
    for (size_t i = 1; i < v1.size(); i++) { p[i] = v1[i] * p[i - 1]; }
    return p;
  }

  template<class T>
  static T sum(const std::vector<T>& v1)
  {
    T s = 0;
    for (const auto& x : v1) { s += x; }
    return s;
  }

  template<class T>
  static T sumProd(const std::vector<T>& v1, const std::vector<T>& v2)
  {
    size_t size;
    if (v1.size() != v2.size())
      throw DimensionException("VectorTools::sumProd", v1.size(), v2.size());
    else
      size = v1.size();

    T x = v2[0] * v1[0];
    for (size_t i = 1; i < size; i++)
    {
      x += v2[i] * v1[i];
    }

    return x;
  }

  template<class T>
  static std::vector<T> cumSum(const std::vector<T>& v1)
  {
    auto s(v1);
    std::partial_sum(s.begin(), s.end(), s.begin());
    return s;
  }

  template<class T>
  static void logNorm(std::vector<T>& v)
  {
    v -= VectorTools::logSumExp(v);
  }

  template<class T>
  static T logSumExp(const std::vector<T>& v1)
  {
    if (v1.size() == 1)
      return v1[0];

    T M = max(v1);
    if (std::isinf(M))
      return M;

    T x = std::accumulate(std::next(v1.begin()), v1.end(), std::exp(v1[0] - M),
          [&M](T y, T z){
        return y + std::exp(z - M);
      });

    return std::log(x) + M;
  }

  template<class T>
  static T logSumExp(const std::vector<T>& v1, const std::vector<T>& v2)
  {
    if (v1.size() != v2.size())
      throw DimensionException("VectorTools::logsumexp", v1.size(), v2.size());

    size_t size = v1.size();

    T M = max(v1);
    if (std::isinf(M))
      throw BadNumberException("VectorTools::logSumExp", M);

    T x = v2[0] * std::exp(v1[0] - M);
    for (size_t i = 1; i < size; i++)
    {
      x += v2[i] * std::exp(v1[i] - M);
    }
    return std::log(x) + M;
  }

  template<class T>
  static T logMeanExp(const std::vector<T>& v1)
  {
    return VectorTools::logSumExp(v1) - std::log(v1.size());
  }


  template<class T>
  static T sumExp(const std::vector<T>& v1)
  {
    if (v1.size() == 0)
      return std::exp(v1[0]);

    T M = max(v1);
    if (std::isinf(M))
      return M < 0 ? 0 : M;

    T x = std::accumulate(std::next(v1.begin()), v1.end(), std::exp(v1[0] - M),
          [&M](T y, T z){
        return y + std::exp(z - M);
      });
    return x * std::exp(M);
  }

  template<class T>
  static T sumExp(const std::vector<T>& v1, const std::vector<T>& v2)
  {
    if (v1.size() != v2.size())
      throw DimensionException("VectorTools::sumExp", v1.size(), v2.size());

    size_t size = v1.size();

    if (size == 1)
      return v2[0] * std::exp(v1[0]);

    T M = max(v1);
    if (std::isinf(M))
      throw BadNumberException("VectorTools::sumExp", M);

    T x = v2[0] * std::exp(v1[0] - M);
    for (size_t i = 1; i < size; i++)
    {
      x += v2[i] * std::exp(v1[i] - M);
    }
    return x * std::exp(M);
  }

  template<class T>
  static std::vector<T> log(const std::vector<T>& v1)
  {
    std::vector<T> v2(v1.size());
    for (size_t i = 0; i < v2.size(); i++)
    {
      v2[i] = log(v1[i]);
    }
    return v2;
  }

  template<class T>
  static double log(const T& x)
  {
    return std::log(x);
  }

  template<class T>
  static std::vector<double> log(const std::vector<T>& v1, double base)
  {
    std::vector<double> v2(v1.size());
    for (size_t i = 0; i < v2.size(); i++) { v2[i] = std::log(v1[i]) / std::log(base); }
    return v2;
  }

  template<class T>
  static double exp(const T& x)
  {
    return std::exp(x);
  }

  template<class T>
  static std::vector<T> exp(const std::vector<T>& v1)
  {
    std::vector<T> v2(v1.size());
    for (size_t i = 0; i < v2.size(); i++) { v2[i] = VectorTools::exp(v1[i]); }
    return v2;
  }

  template<class T>
  static std::vector<double> cos(const std::vector<T>& v1)
  {
    std::vector<double> v2(v1.size());
    for (size_t i = 0; i < v2.size(); i++) { v2[i] = std::cos(v1[i]); }
    return v2;
  }

  template<class T>
  static std::vector<double> sin(const std::vector<T>& v1)
  {
    std::vector<double> v2(v1.size());
    for (size_t i = 0; i < v2.size(); i++) { v2[i] = std::sin(v1[i]); }
    return v2;
  }

  template<class T>
  static std::vector<double> log10(const std::vector<T>& v1)
  {
    std::vector<double> v2(v1.size());
    for (size_t i = 0; i < v1.size(); i++) { v2[i] = std::log10(v1[i]); }
    return v2;
  }

  template<class T>
  static std::vector<T> fact(const std::vector<T>& v1)
  {
    std::vector<T> v2(v1.size());
    for (size_t i = 0; i < v1.size(); i++) { v2[i] = NumTools::fact<T>(v1[i]); }
    return v2;
  }

  template<class T>
  static std::vector<T> sqr(const std::vector<T>& v1)
  {
    std::vector<T> v2(v1.size());
    for (size_t i = 0; i < v1.size(); i++) { v2[i] = NumTools::sqr<T>(v1[i]); }
    return v2;
  }

  template<class T>
  static std::vector<T> pow(const std::vector<T>& v1, T& b)
  {
    std::vector<T> v2(v1.size());
    for (size_t i = 0; i < v1.size(); i++) { v2[i] = std::pow(v1[i], b); }
    return v2;
  }
  template<class T>
  static std::string paste(const std::vector<T>& v, const std::string& delim = " ")
  {
    std::ostringstream out;
    out.precision(12);
    for (size_t i = 0; i < v.size(); i++)
    {
      out << v[i];
      if (i < v.size() - 1)
        out << delim;
    }
    return out.str();
  }

  template<class T>
  static void print(const std::vector<T>& v1, OutputStream& out = * ApplicationTools::message, const std::string& delim = " ")
  {
    for (size_t i = 0; i < v1.size(); i++)
    {
      out << v1[i];
      if (i < v1.size() - 1)
        out << delim;
    }
    out.endLine();
  }

  template<class T>
  static void printRange(const std::vector<T>& v1, OutputStream& out = * ApplicationTools::message, const std::string& delim = ",", const std::string& rangeDelim = "-")
  {
    size_t vs = v1.size();

    for (size_t i = 0; i < vs; i++)
    {
      out << v1[i];
      size_t j = i + 1;

      while (j < vs)
        if (v1[j] == v1[j - 1] + 1)
          j++;
        else
          break;

      if (j > i + 2)
        out << rangeDelim << v1[j - 1];
      i = j - 1;
      if (i < vs - 1)
        out << delim;
    }
  }


  template<class T>
  static void printForR(const std::vector<T>& v1, std::string variableName = "x", std::ostream& out = std::cout)
  {
    out.precision(12);
    out << variableName << "<-c(";
    for (size_t i = 0; i < v1.size(); i++)
    {
      out << v1[i];
      if (i < v1.size() - 1)
        out << ", ";
    }
    out << ")" << std::endl;
  }

  template<class InputType, class OutputType>
  static OutputType scalar(const std::vector<InputType>& v1, const std::vector<InputType>& v2)
  {
    if (v1.size() != v2.size())
    {
      throw DimensionException("VectorTools::scalar", v1.size(), v2.size());
    }
    OutputType result = 0;
    for (size_t i = 0; i < v1.size(); i++)
    {
      result += v1[i] * v2[i];
    }
    return result;
  }
  template<class InputType, class OutputType>
  static OutputType scalar(const std::vector<InputType>& v1, const std::vector<InputType>& v2, const std::vector<InputType>& w)
  {
    if (v1.size() != w.size())
    {
      throw DimensionException("VectorTools::scalar", v1.size(), w.size());
    }
    if (v2.size() != w.size())
    {
      throw DimensionException("VectorTools::scalar", v2.size(), w.size());
    }
    OutputType result = 0;
    for (size_t i = 0; i < v1.size(); i++)
    {
      result += v1[i] * v2[i] * w[i];
    }
    return result;
  }

  template<class T>
  static std::vector<T> kroneckerMult(const std::vector<T>& v1, const std::vector<T>& v2)
  {
    size_t n1 = v1.size();
    size_t n2 = v2.size();
    std::vector<T> v3(n1 * n2);
    for (size_t i = 0; i < n1; i++)
    {
      T v1i = v1[i];
      for (size_t j = 0; j < n2; j++)
      {
        v3[i * n2 + j] = v1i * v2[j];
      }
    }
    return v3;
  }

  template<class InputType, class OutputType>
  static OutputType norm(const std::vector<InputType>& v1)
  {
    OutputType result = 0;
    for (size_t i = 0; i < v1.size(); i++)
    {
      result += v1[i] * v1[i];
    }
    return sqrt(result);
  }

  template<class InputType, class OutputType>
  static OutputType norm(const std::vector<InputType>& v1, const std::vector<InputType>& w)
  {
    if (v1.size() != w.size())
    {
      throw DimensionException("VectorTools::norm", v1.size(), w.size());
    }
    OutputType result = 0;
    for (size_t i = 0; i < v1.size(); i++)
    {
      result += v1[i] * v1[i] * w[i];
    }
    return sqrt(result);
  }

  template<class InputType, class OutputType>
  static OutputType cos(const std::vector<InputType>& v1, const std::vector<InputType>& v2)
  {
    return scalar<InputType, OutputType>(v1, v2)
           / (norm<InputType, OutputType>(v1) * norm<InputType, OutputType>(v2));
  }

  template<class InputType, class OutputType>
  static OutputType cos(const std::vector<InputType>& v1, const std::vector<InputType>& v2, const std::vector<InputType>& w)
  {
    return scalar<InputType, OutputType>(v1, v2, w)
           / (norm<InputType, OutputType>(v1, w) * norm<InputType, OutputType>(v2, w));
  }

  template<class T>
  static T min(const std::vector<T>& v)
  {
    if (v.size() == 0) throw EmptyVectorException<T>("VectorTools::min()", &v);
    T mini = v[0];
    for (size_t i = 1; i < v.size(); i++)
    {
      if (v[i] < mini) mini = v[i];
    }
    return mini;
  }

  template<class T>
  static T max(const std::vector<T>& v)
  {
    if (v.size() == 0) throw EmptyVectorException<T>("VectorTools::max()", &v);
    T maxi = v[0];
    for (size_t i = 1; i < v.size(); i++)
    {
      if (v[i] > maxi) maxi = v[i];
    }
    return maxi;
  }

  template<class T>
  static size_t whichMax(const std::vector<T>& v)
  {
    if (v.size() == 0) throw EmptyVectorException<T>("VectorFuntions::whichMax()", &v);
    T maxi = v[0];
    size_t pos = 0;
    for (size_t i = 1; i < v.size(); i++)
    {
      if (v[i] > maxi)
      {
        maxi = v[i];
        pos = i;
      }
    }
    return pos;
  }

  template<class T>
  static size_t whichMin(const std::vector<T>& v)
  {
    if (v.size() == 0) throw EmptyVectorException<T>("VectorTools::whichMin()", &v);
    T mini = v[0];
    size_t pos = 0;
    for (size_t i = 1; i < v.size(); i++)
    {
      if (v[i] < mini)
      {
        mini = v[i];
        pos = i;
      }
    }
    return pos;
  }

  template<class T>
  static std::vector<size_t> whichMaxAll(const std::vector<T>& v)
  {
    if (v.size() == 0) throw EmptyVectorException<T>("VectorFuntions::whichMaxAll()", &v);
    T maxi = max(v);
    std::vector<size_t> pos;
    for (size_t i = 0; i < v.size(); i++)
    {
      if (v[i] == maxi)
      {
        pos.push_back(i);
      }
    }
    return pos;
  }

  template<class T>
  static std::vector<size_t> whichMinAll(const std::vector<T>& v)
  {
    if (v.size() == 0) throw EmptyVectorException<T>("VectorFuntions::whichMinAll()", &v);
    T mini = min(v);
    std::vector<size_t> pos;
    for (size_t i = 0; i < v.size(); i++)
    {
      if (v[i] == mini)
      {
        pos.push_back(i);
      }
    }
    return pos;
  }


  template<class T>
  static std::vector<T> range(const std::vector<T>& v)
  {
    if (v.size() == 0)
      throw EmptyVectorException<T>("VectorTools::range()", &v);
    std::vector<T> r(2);
    r[0] = r[1] = v[0];
    for (size_t i = 1; i < v.size(); i++)
    {
      if (v[i] < r[0]) r[0] = v[i];
      if (v[i] > r[1]) r[1] = v[i];
    }
    return r;
  }

private:
  template<class T> class order_Cmp_
  {
    const std::vector<T>& values_;

public:
    order_Cmp_(const std::vector<T>& v) : values_(v) {}
    bool operator()(size_t a, size_t b) { return values_[a] < values_[b]; }
  };

public:
  template<class T>
  static std::vector<size_t> order(const std::vector<T>& v)
  {
    if (v.size() == 0)
      throw EmptyVectorException<T>("VectorTools::sort()", &v);
    // Private inner class:
    std::vector<size_t> index(v.size());
    for (size_t i = 0; i < index.size(); ++i)
    {
      index[i] = i;
    }
    sort(index.begin(), index.end(), order_Cmp_<T>(v));
    return index;
  }

  template<class T>
  static std::vector<T> abs(const std::vector<T>& v)
  {
    std::vector<T> vabs(v.size());
    for (size_t i = 0; i < v.size(); i++)
    {
      vabs[i] = std::abs(v[i]);
    }
    return vabs;
  }

  template<class InputType, class OutputType>
  static OutputType mean(const std::vector<InputType>& v1)
  {
    return (OutputType)sum<InputType>(v1) / (OutputType)v1.size();
  }
  template<class InputType, class OutputType>
  static OutputType mean(const std::vector<InputType>& v1, const std::vector<InputType>& w, bool normalizeWeights = true)
  {
    if (normalizeWeights)
    {
      std::vector<InputType> wn = w / sum(w);
      return scalar<InputType, OutputType>(v1, wn);
    }
    else
    {
      return scalar<InputType, OutputType>(v1, w);
    }
  }

  template<class InputType>
  static InputType median(std::vector<InputType>& v1)
  {
    InputType med = 0;
    if (v1.size() == 0) return med;
    if (v1.size() == 1) return v1[0];
    sort(v1.begin(), v1.end());
    size_t i = v1.size() / 2;
    if (v1.size() % 2 == 0)
    {
      // Vector size is pair
      med = double((v1[i - 1] + v1[i]) / 2);
    }
    else
    {
      // Vector size is impair
      med = v1[i];
    }
    return med;
  }

  template<class InputType, class OutputType>
  static std::vector<OutputType> center(const std::vector<InputType>& v1)
  {
    OutputType m = mean<InputType, OutputType>(v1);
    std::vector<OutputType> v(v1.size());
    for (size_t i = 0; i < v1.size(); i++)
    {
      v[i] = (OutputType)v1[i] - m;
    }
    return v;
  }
  template<class InputType, class OutputType>
  static std::vector<OutputType> center(const std::vector<InputType>& v1, const std::vector<InputType>& w, bool normalizeWeights = true)
  {
    OutputType m = mean<InputType, OutputType>(v1, w, normalizeWeights);
    std::vector<OutputType> v(v1.size());
    for (size_t i = 0; i < v1.size(); i++)
    {
      v[i] = (OutputType)v1[i] - m;
    }
    return v;
  }

  template<class InputType, class OutputType>
  static OutputType cov(const std::vector<InputType>& v1, const std::vector<InputType>& v2, bool unbiased = true)
  {
    OutputType n = (OutputType)v1.size();
    OutputType x =  scalar<InputType, OutputType>(
          center<InputType, OutputType>(v1),
          center<InputType, OutputType>(v2)
          ) / n;
    if (unbiased) x = x * n / (n - 1);
    return x;
  }

  template<class InputType, class OutputType>
  static OutputType cov(const std::vector<InputType>& v1, const std::vector<InputType>& v2, const std::vector<InputType>& w, bool unbiased = true, bool normalizeWeights = true)
  {
    if (normalizeWeights)
    {
      std::vector<InputType> wn = w / sum(w);
      OutputType x = scalar<InputType, OutputType>(
            center<InputType, OutputType>(v1, wn, false),
            center<InputType, OutputType>(v2, wn, false),
            wn
            );
      if (unbiased)
      {
        x = x / (1 - sum(sqr<double>(wn)));
      }
      return x;
    }
    else
    {
      OutputType x = scalar<InputType, OutputType>(
            center<InputType, OutputType>(v1, w, false),
            center<InputType, OutputType>(v2, w, false),
            w
            );
      if (unbiased)
      {
        x = x / (1 - sum(sqr(w)));
      }
      return x;
    }
  }
  template<class InputType, class OutputType>
  static OutputType var(const std::vector<InputType>& v1, bool unbiased = true)
  {
    return cov<InputType, OutputType>(v1, v1, unbiased);
  }
  template<class InputType, class OutputType>
  static OutputType var(const std::vector<InputType>& v1, const std::vector<InputType>& w, bool unbiased = true, bool normalizeWeights = true)
  {
    return cov<InputType, OutputType>(v1, v1, w, unbiased, normalizeWeights);
  }

  template<class InputType, class OutputType>
  static OutputType sd(const std::vector<InputType>& v1, bool unbiased = true)
  {
    return sqrt(var<InputType, OutputType>(v1, unbiased));
  }

  template<class InputType, class OutputType>
  static OutputType sd(const std::vector<InputType>& v1, const std::vector<InputType>& w, bool unbiased = true, bool normalizeWeights = true)
  {
    return sqrt(var<InputType, OutputType>(v1, w, unbiased, normalizeWeights));
  }

  template<class InputType, class OutputType>
  static OutputType cor(const std::vector<InputType>& v1, const std::vector<InputType>& v2)
  {
    return cov<InputType, OutputType>(v1, v2)
           / ( sd<InputType, OutputType>(v1) * sd<InputType, OutputType>(v2) );
  }

  template<class InputType, class OutputType>
  static OutputType cor(const std::vector<InputType>& v1, const std::vector<InputType>& v2, const std::vector<InputType>& w, bool normalizeWeights = true)
  {
    if (normalizeWeights)
    {
      std::vector<InputType> wn = w / sum(w);
      return cov<InputType, OutputType>(v1, v2, wn, false, false)
             / ( sd<InputType, OutputType>(v1, wn, false, false) * sd<InputType, OutputType>(v2, wn, false, false) );
    }
    else
    {
      return cov<InputType, OutputType>(v1, v2, w, false, false)
             / ( sd<InputType, OutputType>(v1, w, false, false) * sd<InputType, OutputType>(v2, w, false, false) );
    }
  }

  template<class InputType, class OutputType>
  static OutputType shannon(const std::vector<InputType>& v, double base = 2.7182818)
  {
    OutputType s = 0;
    for (auto x : v)
    {
      if (x > 0)
        s += static_cast<OutputType>(x * std::log(x) / std::log(base));
    }
    return -s;
  }

  template<class InputType, class OutputType>
  static OutputType shannonDiscrete(const std::vector<InputType>& v, double base = 2.7182818)
  {
    std::map<InputType, double> counts;
    for (auto x : v)
    {
      counts[x]++;
    }
    OutputType s = 0;
    double n = static_cast<double>(v.size());
    for (auto& it : counts)
    {
      s += static_cast<OutputType>((it.second / n) * std::log(it.second / n) / std::log(base));
    }
    return -s;
  }

  template<class InputType, class OutputType>
  static OutputType miDiscrete(const std::vector<InputType>& v1, const std::vector<InputType>& v2, double base = 2.7182818)
  {
    if (v1.size() != v2.size())
      throw DimensionException("VectorTools::miDiscrete. The two samples must have the same length.", v2.size(), v1.size());
    std::map<InputType, double> counts1;
    std::map<InputType, double> counts2;
    std::map<InputType, std::map<InputType, double>> counts12;
    for (size_t i = 0; i < v1.size(); i++)
    {
      counts1[v1[i]]++;
      counts2[v2[i]]++;
      counts12[v1[i]][v2[i]]++;
    }
    OutputType s = 0;
    double n = static_cast<double>(v1.size());
    for (auto& it1 : counts12)
    {
      for (auto& it2 : it1.second)
      {
        s += static_cast<OutputType>((it2.second / n) * std::log(it2.second * n / (counts1[it1.first] * counts2[it2.first])) / std::log(base));
      }
    }
    return s;
  }

  template<class InputType, class OutputType>
  static OutputType shannonContinuous(const std::vector<InputType>& v, double base = 2.7182818)
  {
    LinearMatrix<InputType> m(1, v.size());
    for (size_t i = 0; i < v.size(); i++)
    {
      m(0, i) = v[i];
    }
    AdaptiveKernelDensityEstimation kd(m);
    OutputType s = 0;
    std::vector<double> x(1);
    for (auto it : v)
    {
      x[0] = static_cast<double>(it);
      s += static_cast<OutputType>(std::log(kd.kDensity(x)) / std::log(base));
    }
    return -s / static_cast<double>(v.size());
  }

  template<class InputType, class OutputType>
  static OutputType miContinuous(const std::vector<InputType>& v1, const std::vector<InputType>& v2, double base = 2.7182818)
  {
    if (v1.size() != v2.size())
      throw DimensionException("VectorTools::miContinuous. The two samples must have the same length.", v2.size(), v1.size());
    LinearMatrix<InputType> m1(1, v1.size());
    LinearMatrix<InputType> m2(1, v2.size());
    LinearMatrix<InputType> m12(2, v1.size());
    for (size_t i = 0; i < v1.size(); i++)
    {
      m1(0, i) = m12(0, i) = v1[i];
      m2(0, i) = m12(1, i) = v2[i];
    }
    AdaptiveKernelDensityEstimation kd1(m1);
    AdaptiveKernelDensityEstimation kd2(m2);
    AdaptiveKernelDensityEstimation kd12(m12);
    OutputType s = 0;
    std::vector<double> x1(1);
    std::vector<double> x2(1);
    std::vector<double> x12(2);
    for (size_t i = 0; i < v1.size(); i++)
    {
      x1[0] = x12[0] = static_cast<double>(v1[i]);
      x2[0] = x12[1] = static_cast<double>(v2[i]);
      s += static_cast<OutputType>(std::log(kd12.kDensity(x12) / (kd1.kDensity(x1) * kd2.kDensity(x2))) / std::log(base));
    }
    return s / static_cast<double>(v1.size());
  }

  template<class T>
  static bool haveSameElements(const std::vector<T>& v1, const std::vector<T>& v2)
  {
    std::vector<T> u1(v1);
    std::vector<T> u2(v2);
    if (u1.size() != u2.size()) return false;
    std::sort(u1.begin(), u1.end());
    std::sort(u2.begin(), u2.end());
    return u1 == u2;
  }

  template<class T>
  static bool haveSameElements(std::vector<T>& v1, std::vector<T>& v2)
  {
    if (v1.size() != v2.size()) return false;
    std::sort(v1.begin(), v1.end());
    std::sort(v2.begin(), v2.end());
    return v1 == v2;
  }

  template<class T>
  static bool contains(const std::vector<T>& vec, T el)
  {
    for (auto it : vec)
    {
      if (it == el) return true;
    }
    return false;
  }

  template<class T, class U>
  static bool contains(const std::vector<T>& vec, U el)
  {
    for (auto it : vec)
    {
      if (it == (T)el) return true;
    }
    return false;
  }

  template<class T>
  static bool containsAll(std::vector<T>& v1, std::vector<T>& v2)
  {
    std::sort(v1.begin(), v1.end());
    std::sort(v2.begin(), v2.end());
    size_t j = 0;
    for (size_t i = 0; i < v2.size(); i++)
    {
      if (i > 0 && v2[i] == v2[i - 1]) continue;
      while (j < v1.size() - 1 && v1[j] < v2[i]) j++;
      if (v1[j] != v2[i]) return false;
    }
    return true;
  }

  template<class T>
  static std::vector<T> vectorUnion(const std::vector<T>& vec1, const std::vector<T>& vec2)
  {
    std::vector<T> unionEl = vec1;
    for (auto it : vec2)
    {
      if (!contains(unionEl, it))
        unionEl.push_back(it);
    }
    return unionEl;
  }

  template<class T>
  static std::vector<T> vectorUnion(const std::vector< std::vector<T>>& vecElementL)
  {
    std::vector<T> unionEl;
    for (auto it : vecElementL)
    {
      for (auto it2 : it)
      {
        if (!contains(unionEl, it2))
          unionEl.push_back(it2);
      }
    }
    return unionEl;
  }

  template<class T>
  static std::vector<T> vectorIntersection(const std::vector<T>& vec1, const std::vector<T>& vec2)
  {
    std::vector<T> interEl;
    for (auto it : vec1)
    {
      if (contains(vec2, it)) interEl.push_back(it);
    }
    return interEl;
  }

  template<class T, class U>
  static std::vector<T> vectorIntersection(const std::vector<T>& vec1, const std::vector<U>& vec2)
  {
    std::vector<T> interEl;
    for (auto it : vec1)
    {
      if (contains(vec2, it)) interEl.push_back(it);
    }
    return interEl;
  }

  template<class T>
  static std::vector<T> vectorIntersection(const std::vector< std::vector<T>>& vecElementL)
  {
    if (vecElementL.size() == 1) return vecElementL[0];
    std::vector<T> interEl;
    if (vecElementL.size() == 0) return interEl;
    for (auto it : vecElementL[0])
    {
      bool test = true;
      for (size_t j = 1; test && j < vecElementL.size(); j++)
      {
        if (!contains(vecElementL[j], it)) test = false;
      }
      if (test) interEl.push_back(it);
    }
    return interEl;
  }

  template<class T>
  static void append(std::vector<T>& vec1, const std::vector<T>& vec2)
  {
    vec1.insert(vec1.end(), vec2.begin(), vec2.end());
  }

  template<class T>
  static void prepend(std::vector<T>& vec1, const std::vector<T>& vec2)
  {
    vec1.insert(vec1.begin(), vec2.begin(), vec2.end());
  }


  template<class T>
  static std::vector<T> append(const std::vector< std::vector<T>>& vecElementL)
  {
    if (vecElementL.size() == 1) return vecElementL[0];
    std::vector<T> v;
    if (vecElementL.size() == 0) return v;
    for (auto it : vecElementL[0])
    {
      v.push_back(it);
    }
    return v;
  }

  template<class T>
  static void extend(std::vector<T>& vec1, const std::vector<T>& vec2)
  {
    for (auto it : vec2)
    {
      if (!contains(vec1, it))
        vec1.push_back(it);
    }
  }

  template<class T>
  static std::vector<T> rep(const std::vector<T>& vec, size_t n)
  {
    if (n == 1) return vec;
    std::vector<T> v;
    if (n == 0) return v;
    v.resize(vec.size() * n);
    for (size_t i = 0; i < v.size(); i++)
    {
      v[i] = vec[i % vec.size()];
    }
    return v;
  }

  template<class T>
  static void diff(std::vector<T>& v1, std::vector<T>& v2, std::vector<T>& v3)
  {
    if (v2.size() == 0) append(v3, v1);
    std::sort(v1.begin(), v1.end());
    std::sort(v2.begin(), v2.end());
    size_t j = 0;
    for (size_t i = 0; i < v1.size(); i++)
    {
      if (i > 0 && v1[i] == v1[i - 1]) continue;
      while (j < v2.size() - 1 && v2[j] < v1[i]) j++;
      if (v2[j] != v1[i]) v3.push_back(v1[i]);
    }
  }

  static bool test();
};
} // end of namespace bpp.
#endif // BPP_NUMERIC_VECTORTOOLS_H