matfunc.h

Go to the documentation of this file.

#ifndef MATFUNC_H
#define MATFUNC_H

#include <itpp/base/mat.h>
#include <itpp/base/math/log_exp.h>
#include <itpp/base/math/elem_math.h>
#include <itpp/base/algebra/inv.h>
#include <itpp/base/algebra/svd.h>


namespace itpp
{



template<class T>
int length(const Vec<T> &v) { return v.length(); }

template<class T>
int size(const Vec<T> &v) { return v.length(); }

template<class T>
T sum(const Vec<T> &v)
{
  T M = 0;

  for (int i = 0;i < v.length();i++)
    M += v[i];

  return M;
}

template<class T>
Vec<T> sum(const Mat<T> &m, int dim = 1)
{
  it_assert((dim == 1) || (dim == 2), "sum: dimension need to be 1 or 2");
  Vec<T> out;

  if (dim == 1) {
    out.set_size(m.cols(), false);

    for (int i = 0; i < m.cols(); i++)
      out(i) = sum(m.get_col(i));
  }
  else {
    out.set_size(m.rows(), false);

    for (int i = 0; i < m.rows(); i++)
      out(i) = sum(m.get_row(i));
  }

  return out;
}


template<class T>
T sumsum(const Mat<T> &X)
{
  const T * X_data = X._data();
  const int X_datasize = X._datasize();
  T acc = 0;

  for (int i = 0;i < X_datasize;i++)
    acc += X_data[i];

  return acc;
}


template<class T>
T sum_sqr(const Vec<T> &v)
{
  T M = 0;

  for (int i = 0; i < v.length(); i++)
    M += v[i] * v[i];

  return M;
}

template<class T>
Vec<T> sum_sqr(const Mat<T> &m, int dim = 1)
{
  it_assert((dim == 1) || (dim == 2), "sum_sqr: dimension need to be 1 or 2");
  Vec<T> out;

  if (dim == 1) {
    out.set_size(m.cols(), false);

    for (int i = 0; i < m.cols(); i++)
      out(i) = sum_sqr(m.get_col(i));
  }
  else {
    out.set_size(m.rows(), false);

    for (int i = 0; i < m.rows(); i++)
      out(i) = sum_sqr(m.get_row(i));
  }

  return out;
}

template<class T>
Vec<T> cumsum(const Vec<T> &v)
{
  Vec<T> out(v.size());

  out(0) = v(0);
  for (int i = 1; i < v.size(); i++)
    out(i) = out(i - 1) + v(i);

  return out;
}

template<class T>
Mat<T> cumsum(const Mat<T> &m, int dim = 1)
{
  it_assert((dim == 1) || (dim == 2), "cumsum: dimension need to be 1 or 2");
  Mat<T> out(m.rows(), m.cols());

  if (dim == 1) {
    for (int i = 0; i < m.cols(); i++)
      out.set_col(i, cumsum(m.get_col(i)));
  }
  else {
    for (int i = 0; i < m.rows(); i++)
      out.set_row(i, cumsum(m.get_row(i)));
  }

  return out;
}

template<class T>
T prod(const Vec<T> &v)
{
  it_assert(v.size() >= 1, "prod: size of vector should be at least 1");
  T out = v(0);

  for (int i = 1; i < v.size(); i++)
    out *= v(i);

  return out;
}

template<class T>
Vec<T> prod(const Mat<T> &m, int dim = 1)
{
  it_assert((dim == 1) || (dim == 2), "prod: dimension need to be 1 or 2");
  Vec<T> out(m.cols());

  if (dim == 1) {
    it_assert((m.cols() >= 1) && (m.rows() >= 1),
              "prod: number of columns should be at least 1");
    out.set_size(m.cols(), false);

    for (int i = 0; i < m.cols(); i++)
      out(i) = prod(m.get_col(i));
  }
  else {
    it_assert((m.cols() >= 1) && (m.rows() >= 1),
              "prod: number of rows should be at least 1");
    out.set_size(m.rows(), false);

    for (int i = 0; i < m.rows(); i++)
      out(i) = prod(m.get_row(i));
  }
  return out;
}

template<class T>
Vec<T> cross(const Vec<T> &v1, const Vec<T> &v2)
{
  it_assert((v1.size() == 3) && (v2.size() == 3),
            "cross: vectors should be of size 3");

  Vec<T> r(3);

  r(0) = v1(1) * v2(2) - v1(2) * v2(1);
  r(1) = v1(2) * v2(0) - v1(0) * v2(2);
  r(2) = v1(0) * v2(1) - v1(1) * v2(0);

  return r;
}


template<class T>
Vec<T> zero_pad(const Vec<T> &v, int n)
{
  it_assert(n >= v.size(), "zero_pad() cannot shrink the vector!");
  Vec<T> v2(n);
  v2.set_subvector(0, v);
  if (n > v.size())
    v2.set_subvector(v.size(), n - 1, T(0));

  return v2;
}

template<class T>
Vec<T> zero_pad(const Vec<T> &v)
{
  int n = pow2i(levels2bits(v.size()));

  return (n == v.size()) ? v : zero_pad(v, n);
}

template<class T>
Mat<T> zero_pad(const Mat<T> &m, int rows, int cols)
{
  it_assert((rows >= m.rows()) && (cols >= m.cols()),
            "zero_pad() cannot shrink the matrix!");
  Mat<T> m2(rows, cols);
  m2.set_submatrix(0, 0, m);
  if (cols > m.cols()) // Zero
    m2.set_submatrix(0, m.rows() - 1, m.cols(), cols - 1, T(0));
  if (rows > m.rows()) // Zero
    m2.set_submatrix(m.rows(), rows - 1, 0, cols - 1, T(0));

  return m2;
}


template<class T>
T index_zero_pad(const Vec<T> &v, const int index)
{
  if (index >= 0 && index < v.size())
    return v(index);
  else
    return T(0);
}


template<class T>
void transpose(const Mat<T> &m, Mat<T> &out) { out = m.T(); }

template<class T>
Mat<T> transpose(const Mat<T> &m) { return m.T(); }


template<class T>
void hermitian_transpose(const Mat<T> &m, Mat<T> &out) { out = m.H(); }

template<class T>
Mat<T> hermitian_transpose(const Mat<T> &m) { return m.H(); }



template<class Num_T>
bool is_hermitian(const Mat<Num_T>& X)
{

  if (X == X.H())
    return true;
  else
    return false;
}

template<class Num_T>
bool is_unitary(const Mat<Num_T>& X)
{

  if (inv(X) == X.H())
    return true;
  else
    return false;
}


template<class T>
Vec<T> repmat(const Vec<T> &v, int n)
{
  it_assert(n > 0, "repmat(): Wrong repetition parameter");
  int data_length = v.length();
  it_assert(data_length > 0, "repmat(): Input vector can not be empty");
  Vec<T> assembly(data_length * n);
  for (int j = 0; j < n; ++j) {
    assembly.set_subvector(j * data_length, v);
  }
  return assembly;
}


template<class T>
Mat<T> repmat(const Mat<T> &data, int m, int n)
{
  it_assert((m > 0) && (n > 0), "repmat(): Wrong repetition parameters");
  int data_rows = data.rows();
  int data_cols = data.cols();
  it_assert((data_rows > 0) && (data_cols > 0), "repmat(): Input matrix can "
            "not be empty");
  Mat<T> assembly(data_rows*m, data_cols*n);
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) {
      assembly.set_submatrix(i*data_rows, j*data_cols, data);
    }
  }
  return assembly;
}

template<class T> inline
Mat<T> repmat(const Vec<T> &v, int m, int n, bool transpose = false)
{
  return repmat((transpose ? v.T() : Mat<T>(v)), m, n);
}


template<class Num_T>
Mat<Num_T> kron(const Mat<Num_T>& X, const Mat<Num_T>& Y)
{
  Mat<Num_T> result(X.rows() * Y.rows(), X.cols() * Y.cols());

  for (int i = 0; i < X.rows(); i++)
    for (int j = 0; j < X.cols(); j++)
      result.set_submatrix(i * Y.rows(), j * Y.cols(), X(i, j) * Y);

  return result;
}


cmat sqrtm(const cmat& A);

cmat sqrtm(const mat& A);


template<class T>
int rank(const Mat<T> &m, double tol = -1.0)
{
  int rows = m.rows();
  int cols = m.cols();
  if ((rows == 0) || (cols == 0))
    return 0;

  vec sing_val = svd(m);

  if (tol < 0.0) { // Calculate default tolerance
    tol = eps * sing_val(0) * (rows > cols ? rows : cols);
  }

  // Count number of nonzero singular values
  int r = 0;
  while ((r < sing_val.length()) && (sing_val(r) > tol)) {
    r++;
  }

  return r;
}

template<> inline
int rank(const imat &m, double tol)
{
  return rank(to_mat(m), tol);
}

template<> inline
int rank(const smat &m, double tol)
{
  return rank(to_mat(m), tol);
}

template<> inline
int rank(const bmat &, double)
{
  it_error("rank(bmat): Function not implemented for GF(2) algebra");
  return 0;
}




// -------------------- Diagonal matrix functions -------------------------


template<class T>
Mat<T> diag(const Vec<T> &v, const int K = 0)
{
  Mat<T> m(v.size() + std::abs(K), v.size() + std::abs(K));
  m = T(0);
  if (K > 0)
    for (int i = v.size() - 1; i >= 0; i--)
      m(i, i + K) = v(i);
  else
    for (int i = v.size() - 1; i >= 0; i--)
      m(i - K, i) = v(i);

  return m;
}

template<class T>
void diag(const Vec<T> &v, Mat<T> &m)
{
  m.set_size(v.size(), v.size(), false);
  m = T(0);
  for (int i = v.size() - 1; i >= 0; i--)
    m(i, i) = v(i);
}

template<class T>
Vec<T> diag(const Mat<T> &m)
{
  Vec<T> t(std::min(m.rows(), m.cols()));

  for (int i = 0; i < t.size(); i++)
    t(i) = m(i, i);

  return t;
}

template<class T>
Mat<T> bidiag(const Vec<T> &main, const Vec<T> &sup)
{
  it_assert(main.size() == sup.size() + 1, "bidiag()");

  int n = main.size();
  Mat<T> m(n, n);
  m = T(0);
  for (int i = 0; i < n - 1; i++) {
    m(i, i) = main(i);
    m(i, i + 1) = sup(i);
  }
  m(n - 1, n - 1) = main(n - 1);

  return m;
}

template<class T>
void bidiag(const Vec<T> &main, const Vec<T> &sup, Mat<T> &m)
{
  it_assert(main.size() == sup.size() + 1, "bidiag()");

  int n = main.size();
  m.set_size(n, n);
  m = T(0);
  for (int i = 0; i < n - 1; i++) {
    m(i, i) = main(i);
    m(i, i + 1) = sup(i);
  }
  m(n - 1, n - 1) = main(n - 1);
}

template<class T>
void bidiag(const Mat<T> &m, Vec<T> &main, Vec<T> &sup)
{
  it_assert(m.rows() == m.cols(), "bidiag(): Matrix must be square!");

  int n = m.cols();
  main.set_size(n);
  sup.set_size(n - 1);
  for (int i = 0; i < n - 1; i++) {
    main(i) = m(i, i);
    sup(i) = m(i, i + 1);
  }
  main(n - 1) = m(n - 1, n - 1);
}

template<class T>
Mat<T> tridiag(const Vec<T> &main, const Vec<T> &sup, const Vec<T> &sub)
{
  it_assert(main.size() == sup.size() + 1 && main.size() == sub.size() + 1, "bidiag()");

  int n = main.size();
  Mat<T> m(n, n);
  m = T(0);
  for (int i = 0; i < n - 1; i++) {
    m(i, i) = main(i);
    m(i, i + 1) = sup(i);
    m(i + 1, i) = sub(i);
  }
  m(n - 1, n - 1) = main(n - 1);

  return m;
}

template<class T>
void tridiag(const Vec<T> &main, const Vec<T> &sup, const Vec<T> &sub, Mat<T> &m)
{
  it_assert(main.size() == sup.size() + 1 && main.size() == sub.size() + 1, "bidiag()");

  int n = main.size();
  m.set_size(n, n);
  m = T(0);
  for (int i = 0; i < n - 1; i++) {
    m(i, i) = main(i);
    m(i, i + 1) = sup(i);
    m(i + 1, i) = sub(i);
  }
  m(n - 1, n - 1) = main(n - 1);
}

template<class T>
void tridiag(const Mat<T> &m, Vec<T> &main, Vec<T> &sup, Vec<T> &sub)
{
  it_assert(m.rows() == m.cols(), "tridiag(): Matrix must be square!");

  int n = m.cols();
  main.set_size(n);
  sup.set_size(n - 1);
  sub.set_size(n - 1);
  for (int i = 0; i < n - 1; i++) {
    main(i) = m(i, i);
    sup(i) = m(i, i + 1);
    sub(i) = m(i + 1, i);
  }
  main(n - 1) = m(n - 1, n - 1);
}


template<class T>
T trace(const Mat<T> &m)
{
  return sum(diag(m));
}



// ----------------- reshaping vectors and matrices ------------------------


template<class T>
Vec<T> reverse(const Vec<T> &in)
{
  int i, s = in.length();

  Vec<T> out(s);
  for (i = 0;i < s;i++)
    out[i] = in[s-1-i];
  return out;
}

template<class T>
Vec<T> rvectorize(const Mat<T> &m)
{
  int i, j, n = 0, r = m.rows(), c = m.cols();
  Vec<T> v(r * c);

  for (i = 0; i < r; i++)
    for (j = 0; j < c; j++)
      v(n++) = m(i, j);

  return v;
}

template<class T>
Vec<T> cvectorize(const Mat<T> &m)
{
  int i, j, n = 0, r = m.rows(), c = m.cols();
  Vec<T> v(r * c);

  for (j = 0; j < c; j++)
    for (i = 0; i < r; i++)
      v(n++) = m(i, j);

  return v;
}

template<class T>
Mat<T> reshape(const Mat<T> &m, int rows, int cols)
{
  it_assert_debug(m.rows()*m.cols() == rows*cols, "Mat<T>::reshape: Sizes must match");
  Mat<T> temp(rows, cols);
  int i, j, ii = 0, jj = 0;
  for (j = 0; j < m.cols(); j++) {
    for (i = 0; i < m.rows(); i++) {
      temp(ii++, jj) = m(i, j);
      if (ii == rows) {
        jj++;
        ii = 0;
      }
    }
  }
  return temp;
}

template<class T>
Mat<T> reshape(const Vec<T> &v, int rows, int cols)
{
  it_assert_debug(v.size() == rows*cols, "Mat<T>::reshape: Sizes must match");
  Mat<T> temp(rows, cols);
  int i, j, ii = 0;
  for (j = 0; j < cols; j++) {
    for (i = 0; i < rows; i++) {
      temp(i, j) = v(ii++);
    }
  }
  return temp;
}



bool all(const bvec &testvec);
bool any(const bvec &testvec);


// ----------------------------------------------------------------------
// Instantiations
// ----------------------------------------------------------------------

#ifdef HAVE_EXTERN_TEMPLATE

extern template int length(const vec &v);
extern template int length(const cvec &v);
extern template int length(const svec &v);
extern template int length(const ivec &v);
extern template int length(const bvec &v);

extern template double sum(const vec &v);
extern template std::complex<double> sum(const cvec &v);
extern template short sum(const svec &v);
extern template int sum(const ivec &v);
extern template bin sum(const bvec &v);

extern template double sum_sqr(const vec &v);
extern template std::complex<double> sum_sqr(const cvec &v);
extern template short sum_sqr(const svec &v);
extern template int sum_sqr(const ivec &v);
extern template bin sum_sqr(const bvec &v);

extern template vec cumsum(const vec &v);
extern template cvec cumsum(const cvec &v);
extern template svec cumsum(const svec &v);
extern template ivec cumsum(const ivec &v);
extern template bvec cumsum(const bvec &v);

extern template double prod(const vec &v);
extern template std::complex<double> prod(const cvec &v);
extern template short prod(const svec &v);
extern template int prod(const ivec &v);
extern template bin prod(const bvec &v);

extern template vec cross(const vec &v1, const vec &v2);
extern template cvec cross(const cvec &v1, const cvec &v2);
extern template ivec cross(const ivec &v1, const ivec &v2);
extern template svec cross(const svec &v1, const svec &v2);
extern template bvec cross(const bvec &v1, const bvec &v2);

extern template vec reverse(const vec &in);
extern template cvec reverse(const cvec &in);
extern template svec reverse(const svec &in);
extern template ivec reverse(const ivec &in);
extern template bvec reverse(const bvec &in);

extern template vec zero_pad(const vec &v, int n);
extern template cvec zero_pad(const cvec &v, int n);
extern template ivec zero_pad(const ivec &v, int n);
extern template svec zero_pad(const svec &v, int n);
extern template bvec zero_pad(const bvec &v, int n);

extern template vec zero_pad(const vec &v);
extern template cvec zero_pad(const cvec &v);
extern template ivec zero_pad(const ivec &v);
extern template svec zero_pad(const svec &v);
extern template bvec zero_pad(const bvec &v);

extern template mat zero_pad(const mat &, int, int);
extern template cmat zero_pad(const cmat &, int, int);
extern template imat zero_pad(const imat &, int, int);
extern template smat zero_pad(const smat &, int, int);
extern template bmat zero_pad(const bmat &, int, int);

extern template vec sum(const mat &m, int dim);
extern template cvec sum(const cmat &m, int dim);
extern template svec sum(const smat &m, int dim);
extern template ivec sum(const imat &m, int dim);
extern template bvec sum(const bmat &m, int dim);

extern template double sumsum(const mat &X);
extern template std::complex<double> sumsum(const cmat &X);
extern template short sumsum(const smat &X);
extern template int sumsum(const imat &X);
extern template bin sumsum(const bmat &X);

extern template vec sum_sqr(const mat & m, int dim);
extern template cvec sum_sqr(const cmat &m, int dim);
extern template svec sum_sqr(const smat &m, int dim);
extern template ivec sum_sqr(const imat &m, int dim);
extern template bvec sum_sqr(const bmat &m, int dim);

extern template mat cumsum(const mat &m, int dim);
extern template cmat cumsum(const cmat &m, int dim);
extern template smat cumsum(const smat &m, int dim);
extern template imat cumsum(const imat &m, int dim);
extern template bmat cumsum(const bmat &m, int dim);

extern template vec prod(const mat &m, int dim);
extern template cvec prod(const cmat &v, int dim);
extern template svec prod(const smat &m, int dim);
extern template ivec prod(const imat &m, int dim);
extern template bvec prod(const bmat &m, int dim);

extern template vec diag(const mat &in);
extern template cvec diag(const cmat &in);
extern template void diag(const vec &in, mat &m);
extern template void diag(const cvec &in, cmat &m);
extern template mat diag(const vec &v, const int K);
extern template cmat diag(const cvec &v, const int K);

extern template mat bidiag(const vec &, const vec &);
extern template cmat bidiag(const cvec &, const cvec &);
extern template void bidiag(const vec &, const vec &, mat &);
extern template void bidiag(const cvec &, const cvec &, cmat &);
extern template void bidiag(const mat &, vec &, vec &);
extern template void bidiag(const cmat &, cvec &, cvec &);

extern template mat tridiag(const vec &main, const vec &, const vec &);
extern template cmat tridiag(const cvec &main, const cvec &, const cvec &);
extern template void tridiag(const vec &main, const vec &, const vec &, mat &);
extern template void tridiag(const cvec &main, const cvec &, const cvec &, cmat &);
extern template void tridiag(const mat &m, vec &, vec &, vec &);
extern template void tridiag(const cmat &m, cvec &, cvec &, cvec &);

extern template double trace(const mat &in);
extern template std::complex<double> trace(const cmat &in);
extern template short trace(const smat &in);
extern template int trace(const imat &in);
extern template bin trace(const bmat &in);

extern template void transpose(const mat &m, mat &out);
extern template void transpose(const cmat &m, cmat &out);
extern template void transpose(const smat &m, smat &out);
extern template void transpose(const imat &m, imat &out);
extern template void transpose(const bmat &m, bmat &out);

extern template mat transpose(const mat &m);
extern template cmat transpose(const cmat &m);
extern template smat transpose(const smat &m);
extern template imat transpose(const imat &m);
extern template bmat transpose(const bmat &m);

extern template void hermitian_transpose(const mat &m, mat &out);
extern template void hermitian_transpose(const cmat &m, cmat &out);
extern template void hermitian_transpose(const smat &m, smat &out);
extern template void hermitian_transpose(const imat &m, imat &out);
extern template void hermitian_transpose(const bmat &m, bmat &out);

extern template mat hermitian_transpose(const mat &m);
extern template cmat hermitian_transpose(const cmat &m);
extern template smat hermitian_transpose(const smat &m);
extern template imat hermitian_transpose(const imat &m);
extern template bmat hermitian_transpose(const bmat &m);

extern template bool is_hermitian(const mat &X);
extern template bool is_hermitian(const cmat &X);

extern template bool is_unitary(const mat &X);
extern template bool is_unitary(const cmat &X);

extern template vec rvectorize(const mat &m);
extern template cvec rvectorize(const cmat &m);
extern template ivec rvectorize(const imat &m);
extern template svec rvectorize(const smat &m);
extern template bvec rvectorize(const bmat &m);

extern template vec cvectorize(const mat &m);
extern template cvec cvectorize(const cmat &m);
extern template ivec cvectorize(const imat &m);
extern template svec cvectorize(const smat &m);
extern template bvec cvectorize(const bmat &m);

extern template mat reshape(const mat &m, int rows, int cols);
extern template cmat reshape(const cmat &m, int rows, int cols);
extern template imat reshape(const imat &m, int rows, int cols);
extern template smat reshape(const smat &m, int rows, int cols);
extern template bmat reshape(const bmat &m, int rows, int cols);

extern template mat reshape(const vec &m, int rows, int cols);
extern template cmat reshape(const cvec &m, int rows, int cols);
extern template imat reshape(const ivec &m, int rows, int cols);
extern template smat reshape(const svec &m, int rows, int cols);
extern template bmat reshape(const bvec &m, int rows, int cols);

extern template mat kron(const mat &X, const mat &Y);
extern template cmat kron(const cmat &X, const cmat &Y);
extern template imat kron(const imat &X, const imat &Y);
extern template smat kron(const smat &X, const smat &Y);
extern template bmat kron(const bmat &X, const bmat &Y);

extern template vec repmat(const vec &v, int n);
extern template cvec repmat(const cvec &v, int n);
extern template ivec repmat(const ivec &v, int n);
extern template svec repmat(const svec &v, int n);
extern template bvec repmat(const bvec &v, int n);

extern template mat repmat(const vec &v, int m, int n, bool transpose);
extern template cmat repmat(const cvec &v, int m, int n, bool transpose);
extern template imat repmat(const ivec &v, int m, int n, bool transpose);
extern template smat repmat(const svec &v, int m, int n, bool transpose);
extern template bmat repmat(const bvec &v, int m, int n, bool transpose);

extern template mat repmat(const mat &data, int m, int n);
extern template cmat repmat(const cmat &data, int m, int n);
extern template imat repmat(const imat &data, int m, int n);
extern template smat repmat(const smat &data, int m, int n);
extern template bmat repmat(const bmat &data, int m, int n);

#endif // HAVE_EXTERN_TEMPLATE


} // namespace itpp

#endif // #ifndef MATFUNC_H

---