tensor.h

// -*- mode: c++; fill-column: 80; c-basic-offset: 2; indent-tabs-mode: nil -*-
/*
    Copyright (c) 2010 Juan Jose Garcia Ripoll

    Tensor is free software; you can redistribute it and/or modify it
    under the terms of the GNU Library General Public License as published
    by the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Library General Public License for more details.

    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/


#ifndef TENSOR_TENSOR_H
#define TENSOR_TENSOR_H

#include <cassert>
#include <vector>
#include <tensor/numbers.h>
#include <tensor/vector.h>
#include <tensor/gen.h>
#include <tensor/indices.h>
#include <tensor/detail/functional.h>

namespace tensor {

// BASE CLASS
//

/* @{ */
template<typename elt>
class Tensor {
 public:
  typedef elt elt_t;
  typedef elt_t *iterator;
  typedef const elt_t *const_iterator;

  Tensor();

  explicit Tensor(const Indices &new_dims);

  Tensor(const Indices &new_dims, const Tensor<elt_t> &data);

  Tensor(const Vector<elt_t> &data);

  Tensor(const std::vector<elt_t> &data) :
    data_(data.size()), dims_(1)
  {
    dims_.at(0) = data.size();
    std::copy(data.begin(), data.end(), begin());
  }

  Tensor(const Tensor &other);

  template<typename e2> Tensor(const Tensor<e2> &other) :
    data_(other.size()), dims_(other.dimensions())
  {
    std::copy(other.begin(), other.end(), begin());
  }

  template<size_t n> Tensor(const StaticVector<elt_t,n> &t) :
    data_(t), dims_(igen << t.size())
  {}

  template<size_t n> Tensor(const StaticVector<elt_t,n> &t, const Indices &d) :
    data_(t), dims_(d)
  {
    assert(data_.size() == d.total_size());
  }

  explicit Tensor(index length);
  Tensor(index rows, index cols);
  Tensor(index d1, index d2, index d3);
  Tensor(index d1, index d2, index d3, index d4);
  Tensor(index d1, index d2, index d3, index d4, index d5);
  Tensor(index d1, index d2, index d3, index d4, index d5, index d6);

  operator Vector<elt_t>() const { return data_; }

  const Tensor &operator=(const Tensor<elt_t> &other);

  index size() const { return data_.size(); }
  bool is_empty() const { return size() == 0; }

  int rank() const { return dims_.size(); }
  const Indices &dimensions() const { return dims_; }
  index dimension(int which) const;
  void get_dimensions(index *length) const;
  void get_dimensions(index *rows, index *cols) const;
  void get_dimensions(index *d0, index *d1, index *d2) const;
  void get_dimensions(index *d0, index *d1, index *d2, index *d3) const;
  void get_dimensions(index *d0, index *d1, index *d2, index *d3, index *d4) const;
  void get_dimensions(index *d0, index *d1, index *d2, index *d3, index *d4, index *d5) const;
  index columns() const { return dimension(1); }
  index rows() const { return dimension(0); }

  void reshape(const Indices &new_dims);

  const elt_t &operator[](index i) const;
  const elt_t &operator()(index i) const;
  const elt_t &operator()(index row, index col) const;
  const elt_t &operator()(index d0, index d1, index d2) const;
  const elt_t &operator()(index d0, index d1, index d2, index d3) const;
  const elt_t &operator()(index d0, index d1, index d2, index d3, index d4) const;
  const elt_t &operator()(index d0, index d1, index d2, index d3, index d4, index d5w) const;

  elt_t &at_seq(index i);
  elt_t &at(index i);
  elt_t &at(index row, index col);
  elt_t &at(index d1, index d2, index d3);
  elt_t &at(index d1, index d2, index d3, index d4);
  elt_t &at(index d1, index d2, index d3, index d4, index d5);
  elt_t &at(index d1, index d2, index d3, index d4, index d5, index d6);

  void fill_with(const elt_t &e);
  void fill_with_zeros() { fill_with(number_zero<elt_t>()); }
  void randomize();

  static const Tensor<elt_t> random(index length);
  static const Tensor<elt_t> random(index rows, index cols);
  static const Tensor<elt_t> random(index d1, index d2, index d3);
  static const Tensor<elt_t> random(index d1, index d2, index d3, index d4);
  static const Tensor<elt_t> random(index d1, index d2, index d3, index d4, index d5);
  static const Tensor<elt_t> random(index d1, index d2, index d3, index d4, index d5, index d6);
  static const Tensor<elt_t> random(const Indices &dimensions);

  //
  // Tensor slicing
  //
  class view;
  const view operator()(PRange r) const;
  const view operator()(PRange r1, PRange r2) const;
  const view operator()(PRange r1, PRange r2, PRange r3) const;
  const view operator()(PRange r1, PRange r2, PRange r3, PRange r4) const;
  const view operator()(PRange r1, PRange r2, PRange r3, PRange r4, PRange r5) const;
  const view operator()(PRange r1, PRange r2, PRange r3, PRange r4, PRange r5, PRange r6) const;

  class mutable_view;
  mutable_view at(PRange r);
  mutable_view at(PRange r1, PRange r2);
  mutable_view at(PRange r1, PRange r2, PRange r3);
  mutable_view at(PRange r1, PRange r2, PRange r3, PRange r4);
  mutable_view at(PRange r1, PRange r2, PRange r3, PRange r4, PRange r5);
  mutable_view at(PRange r1, PRange r2, PRange r3, PRange r4, PRange r5, PRange r6);

  //
  // Matrix operations
  //
  static const Tensor<elt_t> eye(index rows) { return eye(rows, rows); }
  static const Tensor<elt_t> eye(index rows, index cols);
  static const Tensor<elt_t> zeros(index rows) { return zeros(rows, rows); }
  static const Tensor<elt_t> zeros(index rows, index cols);
  static const Tensor<elt_t> zeros(index d1, index d2, index d3);
  static const Tensor<elt_t> zeros(index d1, index d2, index d3, index d4);
  static const Tensor<elt_t> zeros(index d1, index d2, index d3, index d4, index d5);
  static const Tensor<elt_t> zeros(index d1, index d2, index d3, index d4, index d5, index d6);
  static const Tensor<elt_t> zeros(const Indices &dimensions);
  static const Tensor<elt_t> ones(index rows) { return ones(rows, rows); }
  static const Tensor<elt_t> ones(index rows, index cols);
  static const Tensor<elt_t> ones(const Indices &dimensions);

  iterator begin() { return data_.begin(); }
  const_iterator begin() const { return data_.begin_const(); }
  const_iterator begin_const() const { return data_.begin_const(); }
  const_iterator end_const() const { return data_.end_const(); }
  const_iterator end() const { return data_.end_const(); }
  iterator end() { return data_.end(); }

  // Only for testing purposes
  int ref_count() const { return data_.ref_count(); }

  const Tensor<elt_t> diag(int which = 0, int ndx1 = 0, int ndx2 = -1) {
    return take_diag(*this, which, ndx1, ndx2);
  }

 private:
  Vector<elt_t> data_;
  Indices dims_;
};

// ALGEBRA
//
//
// Unary operations
//
template<typename t>
Tensor<t> operator-(const Tensor<t> &);

//
// Binary operations
//
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator+(const Tensor<t1> &a, const Tensor<t2> &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator-(const Tensor<t1> &a, const Tensor<t2> &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator*(const Tensor<t1> &a, const Tensor<t2> &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator/(const Tensor<t1> &a, const Tensor<t2> &b);

template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator+(const Tensor<t1> &a, const t2 &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator-(const Tensor<t1> &a, const t2 &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator*(const Tensor<t1> &a, const t2 &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator/(const Tensor<t1> &a, const t2 &b);

template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator+(const t1 &a, const Tensor<t2> &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator-(const t1 &a, const Tensor<t2> &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator*(const t1 &a, const Tensor<t2> &b);
template<typename t1, typename t2>
Tensor<typename Binop<t1,t2>::type> operator/(const t1 &a, const Tensor<t2> &b);

  template<typename elt_t>
  const Tensor<elt_t> kron(const Tensor<elt_t> &a, const Tensor<elt_t> &b);
  template<typename elt_t>
  const Tensor<elt_t> kron2(const Tensor<elt_t> &a, const Tensor<elt_t> &b);
  template<typename elt_t>
  const Tensor<elt_t> kron2_sum(const Tensor<elt_t> &a, const Tensor<elt_t> &b);

} // namespace tensor

// IMPLEMENTATIONS
//
#ifdef TENSOR_LOAD_IMPL
#include <tensor/detail/tensor_base.hpp>
#include <tensor/detail/tensor_matrix.hpp>
#include <tensor/detail/tensor_kron.hpp>
#endif
#include <tensor/detail/tensor_slice.hpp>
#include <tensor/detail/tensor_ops.hpp>

// EXPLICIT INSTANTIATIONS
//

namespace tensor {

  extern template class Tensor<double>;
#ifdef DOXYGEN_ONLY
  struct RTensor : public Tensor<double> {}
#else
  typedef Tensor<double> RTensor;
#endif

  const RTensor change_dimension(const RTensor &U, int dimension, index new_size);

  double sum(const RTensor &r);
  double mean(const RTensor &r);
  const RTensor mean(const RTensor &r, int which);

  double norm0(const RTensor &r);
  double scprod(const RTensor &a, const RTensor &b);
  double norm2(const RTensor &r);
  double matrix_norminf(const RTensor &r);

  RTensor abs(const RTensor &t);
  RTensor cos(const RTensor &t);
  RTensor sin(const RTensor &t);
  RTensor tan(const RTensor &t);
  RTensor cosh(const RTensor &t);
  RTensor sinh(const RTensor &t);
  RTensor tanh(const RTensor &t);
  RTensor exp(const RTensor &t);
  RTensor sqrt(const RTensor &t);

  RTensor round(const RTensor &t);

  const RTensor diag(const RTensor &d, int which, int rows, int cols);
  const RTensor diag(const RTensor &d, int which = 0);
  const RTensor take_diag(const RTensor &d, int which = 0, int ndx1 = 0, int ndx2 = -1);
  double trace(const RTensor &d);
  const RTensor trace(const RTensor &A, index i1, index i2);

  const RTensor reshape(const RTensor &t, const Indices &new_dims);
  const RTensor reshape(const RTensor &t, index length);
  const RTensor reshape(const RTensor &t, index rows, index columns);
  const RTensor reshape(const RTensor &t, index d1, index d2, index d3);
  const RTensor reshape(const RTensor &t, index d1, index d2, index d3, index d4);
  const RTensor reshape(const RTensor &t, index d1, index d2, index d3, index d4, index d5);
  const RTensor reshape(const RTensor &t, index d1, index d2, index d3, index d4, index d5, index d6);

  const RTensor squeeze(const RTensor &t);
  const RTensor permute(const RTensor &a, index ndx1 = 0, index ndx2 = -1);
  const RTensor transpose(const RTensor &a);
  inline const RTensor adjoint(const RTensor &a) { return transpose(a); }

  const RTensor fold(const RTensor &a, int ndx1, const RTensor &b, int ndx2);
  const RTensor foldc(const RTensor &a, int ndx1, const RTensor &b, int ndx2);
  const RTensor foldin(const RTensor &a, int ndx1, const RTensor &b, int ndx2);
  const RTensor mmult(const RTensor &a, const RTensor &b);

  void fold_into(RTensor &output, const RTensor &a, int ndx1, const RTensor &b, int ndx2);
  void foldin_into(RTensor &output, const RTensor &a, int ndx1, const RTensor &b, int ndx2);
  void mmult_into(RTensor &output, const RTensor &a, const RTensor &b);

  bool all_equal(const RTensor &a, const RTensor &b);
  bool all_equal(const RTensor &a, double b);
  inline bool all_equal(double b, const RTensor &a) { return all_equal(a, b); }
  template<typename t1, typename t2>
  inline bool some_unequal(const t1 &a, const t2 &b) { return !all_equal(a,b); }

  const Booleans operator==(const RTensor &a, const RTensor &b);
  const Booleans operator<(const RTensor &a, const RTensor &b);
  const Booleans operator>(const RTensor &a, const RTensor &b);
  const Booleans operator<=(const RTensor &a, const RTensor &b);
  const Booleans operator>=(const RTensor &a, const RTensor &b);
  const Booleans operator!=(const RTensor &a, const RTensor &b);

  const Booleans operator==(const RTensor &a, double b);
  const Booleans operator<(const RTensor &a, double b);
  const Booleans operator>(const RTensor &a, double b);
  const Booleans operator<=(const RTensor &a, double b);
  const Booleans operator>=(const RTensor &a, double b);
  const Booleans operator!=(const RTensor &a, double b);

  inline const Booleans operator==(double a, const RTensor &b) { return b == a; }
  inline const Booleans operator<(double a, const RTensor &b) { return b >= a; }
  inline const Booleans operator>(double a, const RTensor &b) { return b <= a; }
  inline const Booleans operator<=(double a, const RTensor &b) { return b > a; }
  inline const Booleans operator>=(double a, const RTensor &b) {  return b < a; }
  inline const Booleans operator!=(double a, const RTensor &b) { return b != a; }

  extern template class Tensor<cdouble>;
#ifdef DOXYGEN_ONLY
  struct CTensor : public Tensor<cdouble> {}
#else
  typedef Tensor<cdouble> CTensor;
#endif

  const CTensor change_dimension(const CTensor &U, int dimension, index new_size);

  cdouble sum(const CTensor &r);
  cdouble mean(const CTensor &r);
  const CTensor mean(const CTensor &r, int ndx);

  double norm0(const CTensor &r);
  cdouble scprod(const CTensor &a, const CTensor &b);
  double norm2(const CTensor &r);
  double matrix_norminf(const CTensor &r);

  inline const RTensor real(const RTensor &r) { return r; }
  const RTensor imag(const RTensor &r);
  const RTensor real(const CTensor &r);
  const RTensor imag(const CTensor &r);

  const CTensor to_complex(const RTensor &r);
  inline const CTensor to_complex(const CTensor &r) { return r; }
  const CTensor to_complex(const RTensor &r, const RTensor &i);

  inline const RTensor conj(const RTensor &r) { return r; }
  const CTensor conj(const CTensor &c);

  RTensor abs(const CTensor &t);
  CTensor cos(const CTensor &t);
  CTensor sin(const CTensor &t);
  CTensor tan(const CTensor &t);
  CTensor cosh(const CTensor &t);
  CTensor sinh(const CTensor &t);
  CTensor tanh(const CTensor &t);
  CTensor exp(const CTensor &t);
  CTensor sqrt(const CTensor &t);

  const CTensor diag(const CTensor &d, int which, int rows, int cols);
  const CTensor diag(const CTensor &d, int which = 0);
  const CTensor take_diag(const CTensor &d, int which = 0, int ndx1 = 0, int ndx2 = -1);
  cdouble trace(const CTensor &d);
  const CTensor trace(const CTensor &A, index i1, index i2);

  const CTensor reshape(const CTensor &t, const Indices &new_dims);
  const CTensor reshape(const CTensor &t, index length);
  const CTensor reshape(const CTensor &t, index rows, index columns);
  const CTensor reshape(const CTensor &t, index d1, index d2, index d3);
  const CTensor reshape(const CTensor &t, index d1, index d2, index d3, index d4);
  const CTensor reshape(const CTensor &t, index d1, index d2, index d3, index d4, index d5);
  const CTensor reshape(const CTensor &t, index d1, index d2, index d3, index d4, index d5, index d6);

  const CTensor squeeze(const CTensor &t);
  const CTensor permute(const CTensor &a, index ndx1 = 0, index ndx2 = -1);
  const CTensor transpose(const CTensor &a);
  const CTensor adjoint(const CTensor &a);

  const CTensor fold(const CTensor &a, int ndx1, const CTensor &b, int ndx2);
  const CTensor fold(const RTensor &a, int ndx1, const CTensor &b, int ndx2);
  const CTensor fold(const CTensor &a, int ndx1, const RTensor &b, int ndx2);

  const CTensor foldc(const CTensor &a, int ndx1, const CTensor &b, int ndx2);
  const CTensor foldc(const RTensor &a, int ndx1, const CTensor &b, int ndx2);
  const CTensor foldc(const CTensor &a, int ndx1, const RTensor &b, int ndx2);

  const CTensor mmult(const CTensor &a, const CTensor &b);
  const CTensor mmult(const RTensor &a, const CTensor &b);
  const CTensor mmult(const CTensor &a, const RTensor &b);

  const RTensor scale(const RTensor &t, int ndx1, const RTensor &v);
  const CTensor scale(const CTensor &t, int ndx1, const CTensor &v);
  const CTensor scale(const CTensor &t, int ndx1, const RTensor &v);
  void scale_inplace(RTensor &t, int ndx1, const RTensor &v);
  void scale_inplace(CTensor &t, int ndx1, const CTensor &v);
  void scale_inplace(CTensor &t, int ndx1, const RTensor &v);

  const CTensor foldin(const CTensor &a, int ndx1, const CTensor &b, int ndx2);

  const RTensor linspace(double min, double max, index n = 100);
  const RTensor linspace(const RTensor &min, const RTensor &max, index n = 100);
  const CTensor linspace(cdouble min, cdouble max, index n = 100);
  const CTensor linspace(const CTensor &min, const CTensor &max, index n = 100);

  const RTensor sort(const RTensor &v, bool reverse = false);
  const Indices sort_indices(const RTensor &v, bool reverse = false);

  const CTensor sort(const CTensor &v, bool reverse = false);
  const Indices sort_indices(const CTensor &v, bool reverse = false);

  bool all_equal(const CTensor &a, const CTensor &b);
  bool all_equal(const CTensor &a, const cdouble &b);
  inline bool all_equal(cdouble b, const CTensor &a) { return all_equal(a, b); }

  const Booleans operator==(const CTensor &a, const CTensor &b);
  const Booleans operator!=(const CTensor &a, const CTensor &b);
  const Booleans operator==(const CTensor &a, cdouble b);
  const Booleans operator!=(const CTensor &a, cdouble b);
  inline const Booleans operator==(cdouble a, const CTensor &b) { return b == a; }
  inline const Booleans operator!=(cdouble a, const CTensor &b) { return b != a; }

} // namespace tensor

/* @} */

#endif // !TENSOR_H