VectorTools.h

Go to the documentation of this file.

//
// File: VectorTools.h
// Authors:
//   Julien Dutheil
// Created: 2003-03-14 14:16:32
//
 
/*
  Copyright or © or Copr. Bio++ Development Team, (November 17, 2004)
  
  This software is a computer program whose purpose is to provide classes
  for numerical calculus.
  
  This software is governed by the CeCILL license under French law and
  abiding by the rules of distribution of free software. You can use,
  modify and/ or redistribute the software under the terms of the CeCILL
  license as circulated by CEA, CNRS and INRIA at the following URL
  "http://www.cecill.info".
  
  As a counterpart to the access to the source code and rights to copy,
  modify and redistribute granted by the license, users are provided only
  with a limited warranty and the software's author, the holder of the
  economic rights, and the successive licensors have only limited
  liability.
  
  In this respect, the user's attention is drawn to the risks associated
  with loading, using, modifying and/or developing or reproducing the
  software by the user in light of its specific status of free software,
  that may mean that it is complicated to manipulate, and that also
  therefore means that it is reserved for developers and experienced
  professionals having in-depth computer knowledge. Users are therefore
  encouraged to load and test the software's suitability as regards their
  requirements in conditions enabling the security of their systems and/or
  data to be ensured and, more generally, to use and operate it in the
  same conditions as regards security.
  
  The fact that you are presently reading this means that you have had
  knowledge of the CeCILL license and that you accept its terms.
*/
 
#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<int>& 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