Tensor.hpp

Go to the documentation of this file.

#ifndef CYTNX_TENSOR_H_
#define CYTNX_TENSOR_H_
 
#include "Type.hpp"
#include "cytnx_error.hpp"
#include "Device.hpp"
#include "intrusive_ptr_base.hpp"
#include <iostream>
#include <fstream>
#include "utils/dynamic_arg_resolver.hpp"
#include "Accessor.hpp"
#include <type_traits>
#include <utility>
#include <vector>
#include <initializer_list>
#include <string>
 
#include "backend/Scalar.hpp"
#include "backend/Storage.hpp"
#include "backend/Tensor_impl.hpp"
 
namespace cytnx {
 
  class Tensor;
 
  // [Note] The free `Tensor <op> T` operators declared in linalg.hpp are forward-declared below
  // the class definition rather than here, because their constraint names Tensor::Tproxy, which
  // is not available until Tensor is complete. Declaring them unconstrained here would defeat the
  // constraint entirely: a differently-constrained declaration is a *distinct* template, so the
  // unconstrained candidate would stay viable for every non-scalar T (#1003, Ian's review).
 
  class Tensor {
   private:
   public:
    // this is a proxy class to allow get/set element using [] as python!
    struct Tproxy {
      boost::intrusive_ptr<Tensor_impl> _insimpl;
      std::vector<cytnx::Accessor> _accs;
      Tproxy(boost::intrusive_ptr<Tensor_impl> _ptr, const std::vector<cytnx::Accessor> &accs)
          : _insimpl(std::move(_ptr)), _accs(accs) {}
 
      // when used to set elems:
      const Tensor &operator=(const Tensor &rhs) {
        this->_insimpl->set(_accs, rhs._impl);
        return rhs;
      }
 
      template <class T>
      const T &operator=(const T &rc) {
        this->_insimpl->set(_accs, rc);
        return rc;
      }
      const Tproxy &operator=(const Tproxy &rc) {
        Tensor tmp = Tensor(rc);
        this->_insimpl->set(_accs, tmp._impl);
        return rc;
      }
 
      template <class T>
      Tensor operator+=(const T &rc) {
        Tensor self;
        self._impl = _insimpl->get(_accs);
        self += rc;
        _insimpl->set(_accs, self._impl);
        self._impl = this->_insimpl;
        return self;
      }
      Tensor operator+=(const Tproxy &rc);
 
      template <class T>
      Tensor operator-=(const T &rc) {
        Tensor self;
        self._impl = _insimpl->get(_accs);
        self -= rc;
        _insimpl->set(_accs, self._impl);
        self._impl = this->_insimpl;
        return self;
      }
      Tensor operator-=(const Tproxy &rc);
 
      template <class T>
      Tensor operator/=(const T &rc) {
        Tensor self;
        self._impl = _insimpl->get(_accs);
        self /= rc;
        _insimpl->set(_accs, self._impl);
        self._impl = this->_insimpl;
        return self;
      }
      Tensor operator/=(const Tproxy &rc);
 
      template <class T>
      Tensor operator*=(const T &rc) {
        Tensor self;
        self._impl = _insimpl->get(_accs);
        self *= rc;
        _insimpl->set(_accs, self._impl);
        self._impl = this->_insimpl;
        return self;
      }
      Tensor operator*=(const Tproxy &rc);
 
      // alias to resolve conflict with op ovld for rc=Tensor
      /*
      template<class T>
      Tensor _operatorADD(const T &rc) const{
          Tensor out;
          out._impl = _insimpl->get(_accs);
          return out.Add(rc);
      }
      */
      Tensor operator+(const cytnx_complex128 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_complex64 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_double &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_float &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_uint64 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_int64 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_uint32 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_int32 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_uint16 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_int16 &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const cytnx_bool &rc) const;  //{return this->_operatorADD(rc);};
      Tensor operator+(const Tproxy &rc) const;
 
      /*
      template<class T>
      Tensor _operatorSUB(const T &rc) const{
          Tensor out;
          out._impl = _insimpl->get(_accs);
          return out.Sub(rc);
      }
      */
      Tensor operator-(const cytnx_complex128 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_complex64 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_double &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_float &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_uint64 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_int64 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_uint32 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_int32 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_uint16 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_int16 &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const cytnx_bool &rc) const;  //{return this->_operatorSUB(rc);};
      Tensor operator-(const Tproxy &rc) const;
 
      Tensor operator-() const;
 
      /*
      template<class T>
      Tensor _operatorMUL(const T &rc) const{
          Tensor out;
          out._impl = _insimpl->get(_accs);
          return out.Mul(rc);
      }
      */
      Tensor operator*(const cytnx_complex128 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_complex64 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_double &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_float &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_uint64 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_int64 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_uint32 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_int32 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_uint16 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_int16 &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const cytnx_bool &rc) const;  //{return this->_operatorMUL(rc);};
      Tensor operator*(const Tproxy &rc) const;
 
      /*
      template<class T>
      Tensor _operatorDIV(const T &rc) const{
          Tensor out;
          out._impl = _insimpl->get(_accs);
          return out.Div(rc);
      }
      */
      Tensor operator/(const cytnx_complex128 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_complex64 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_double &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_float &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_uint64 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_int64 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_uint32 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_int32 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_uint16 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_int16 &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const cytnx_bool &rc) const;  //{return this->_operatorDIV(rc);};
      Tensor operator/(const Tproxy &rc) const;
 
      template <class T>
      T item() const {
        Tensor out;
        out._impl = _insimpl->get(_accs);
        return out.item<T>();
      }
 
      Scalar::Sproxy item() const {
        Tensor out;
        out._impl = _insimpl->get(_accs);
        return out.item();
      }
 
      // when used to get elems:
      operator Tensor() const {
        Tensor out;
        out._impl = _insimpl->get(_accs);
        return out;
      }
 
      Storage storage() const {
        Tensor out;
        out._impl = _insimpl->get(_accs);
        return out.storage();
      }
 
    };  // proxy class of Tensor.
 
 
    // these two are using the python way!
    //----------------------------------------
    template <class... Ts>
    Tproxy operator()(const std::string &e1, const Ts &...elems) {
      std::vector<cytnx::Accessor> tmp = Indices_resolver(e1, elems...);
      return (*this)[tmp];
    }
    template <class... Ts>
    Tproxy operator()(const cytnx_int64 &e1, const Ts &...elems) {
      std::vector<cytnx::Accessor> tmp = Indices_resolver(e1, elems...);
      return (*this)[tmp];
    }
    template <class... Ts>
    Tproxy operator()(const cytnx::Accessor &e1, const Ts &...elems) {
      std::vector<cytnx::Accessor> tmp = Indices_resolver(e1, elems...);
      return (*this)[tmp];
    }
    template <class... Ts>
    const Tproxy operator()(const std::string &e1, const Ts &...elems) const {
      std::vector<cytnx::Accessor> tmp = Indices_resolver(e1, elems...);
      return (*this)[tmp];
    }
    template <class... Ts>
    const Tproxy operator()(const cytnx_int64 &e1, const Ts &...elems) const {
      std::vector<cytnx::Accessor> tmp = Indices_resolver(e1, elems...);
      return (*this)[tmp];
    }
    template <class... Ts>
    const Tproxy operator()(const cytnx::Accessor &e1, const Ts &...elems) const {
      std::vector<cytnx::Accessor> tmp = Indices_resolver(e1, elems...);
      return (*this)[tmp];
    }
 
    //-----------------------------------------
 
    Tproxy operator[](const std::initializer_list<cytnx::Accessor> &accs) {
      std::vector<cytnx::Accessor> tmp = accs;
      return (*this)[tmp];
    }
    Tproxy operator[](const std::vector<cytnx::Accessor> &accs) {
      return Tproxy(this->_impl, accs);
    }
 
    const Tproxy operator[](const std::vector<cytnx::Accessor> &accs) const {
      return Tproxy(this->_impl, accs);
    }
    const Tproxy operator[](const std::initializer_list<cytnx::Accessor> &accs) const {
      std::vector<cytnx::Accessor> tmp = accs;
      return (*this)[tmp];
    }
 
    Tproxy operator[](const std::initializer_list<cytnx_int64> &accs) {
      std::vector<cytnx_int64> tmp = accs;
      return (*this)[tmp];
    }
    Tproxy operator[](const std::vector<cytnx_int64> &accs) {
      std::vector<cytnx::Accessor> acc_in;
      for (int i = 0; i < accs.size(); i++) {
        acc_in.push_back(cytnx::Accessor(accs[i]));
      }
      return Tproxy(this->_impl, acc_in);
    }
    const Tproxy operator[](const std::initializer_list<cytnx_int64> &accs) const {
      std::vector<cytnx_int64> tmp = accs;
      return (*this)[tmp];
    }
    const Tproxy operator[](const std::vector<cytnx_uint64> &accs) const {
      std::vector<cytnx::Accessor> acc_in;
      for (int i = 0; i < accs.size(); i++) {
        acc_in.push_back(cytnx::Accessor(accs[i]));
      }
      return Tproxy(this->_impl, acc_in);
    }
    const Tproxy operator[](const std::vector<cytnx_int64> &accs) const {
      std::vector<cytnx::Accessor> acc_in;
      for (int i = 0; i < accs.size(); i++) {
        acc_in.push_back(cytnx::Accessor(accs[i]));
      }
      return Tproxy(this->_impl, acc_in);
    }
    //-------------------------------------------
 
    void _Save(std::fstream &f) const;
    void _Load(std::fstream &f);
 
 
    void Save(const std::string &fname) const;
    void Save(const char *fname) const;
 
    void Tofile(const std::string &fname) const;
 
    void Tofile(const char *fname) const;
 
    void Tofile(std::fstream &f) const;
 
    static Tensor Load(const std::string &fname);
    static Tensor Load(const char *fname);
 
    static Tensor Fromfile(const std::string &fname, const unsigned int &dtype,
                           const cytnx_int64 &count = -1);
    static Tensor Fromfile(const char *fname, const unsigned int &dtype,
                           const cytnx_int64 &count = -1);
 
    // static Tensor Frombinary(const std::string &fname);
 
    boost::intrusive_ptr<Tensor_impl> _impl;
    Tensor() : _impl(new Tensor_impl()){};
    Tensor(const Tensor &rhs) { _impl = rhs._impl; }
 
    /*
    template<class Tp>
    Tensor(const std::initializer_list<Tp> &rhs){
        Storage stmp = std::vector<Tp>(rhs);
        boost::intrusive_ptr<Tensor_impl> tmp(new Tensor_impl());
        tmp->Init(stmp);
        this->_impl = tmp;
    }
    */
 
    Tensor &operator=(const Tensor &rhs) {
      _impl = rhs._impl;
      return *this;
    }
 
    void operator=(const Tproxy &rhsp) {  // this is used to handle proxy assignment
      this->_impl = rhsp._insimpl->get(rhsp._accs);
    }
 
    // default device==Device.cpu (-1)
    void Init(const std::vector<cytnx_uint64> &shape, unsigned int dtype = Type.Double,
              int device = -1, bool init_zero = true) {
      boost::intrusive_ptr<Tensor_impl> tmp(new Tensor_impl());
      this->_impl = tmp;
      this->_impl->Init(shape, dtype, device, init_zero);
    }
    void Init(std::initializer_list<cytnx_uint64> shape, unsigned int dtype = Type.Double,
              int device = -1, bool init_zero = true) {
      this->Init(std::vector<cytnx_uint64>(shape), dtype, device, init_zero);
    }
    // void Init(const Storage& storage) {
    //   boost::intrusive_ptr<Tensor_impl> tmp(new Tensor_impl());
    //   this->_impl = tmp;
    //   this->_impl->Init(storage);
    // }
    // void Init(const Storage& storage, const std::vector<cytnx_uint64> &shape,
    //   const unsigned int &dtype = Type.Double, const int &device = -1) {
    //   boost::intrusive_ptr<Tensor_impl> tmp(new Tensor_impl());
    //   this->_impl = tmp;
    //   this->_impl->Init(storage, shape, dtype, device);
    // }
 
    Tensor(const std::vector<cytnx_uint64> &shape, unsigned int dtype = Type.Double,
           int device = -1, bool init_zero = true)
        : _impl(new Tensor_impl()) {
      this->Init(shape, dtype, device, init_zero);
    }
    Tensor(std::initializer_list<cytnx_uint64> shape, unsigned int dtype = Type.Double,
           int device = -1, bool init_zero = true)
        : _impl(new Tensor_impl()) {
      this->Init(shape, dtype, device, init_zero);
    }
    // Tensor(const Storage& storage)
    //     : _impl(new Tensor_impl()) {
    //   this->Init(storage);
    // }
    // Tensor(const Storage& storage, const std::vector<cytnx_uint64> &shape,
    //   const unsigned int &dtype = Type.Double, const int &device = -1)
    //     : _impl(new Tensor_impl()) {
    //   this->Init(storage, shape, dtype, device);
    // }
 
    // This mechanism is to remove the 'void' type from Type_list. Taking advantage of it
    // appearing first ...
 
    struct internal {
      template <typename Variant>
      struct exclude_first;
 
      template <typename First, typename... Rest>
      struct exclude_first<std::variant<First, Rest...>> {
        using type = std::variant<Rest...>;
      };
    };  // internal
 
    // std::variant of pointers to Type_list, without void ....
    using pointer_types =
      make_variant_from_transform_t<typename internal::exclude_first<Type_list>::type,
                                    std::add_pointer>;
 
    // convert this->_impl->_storage._impl->Mem to a typed variant of pointers, excluding void*
    pointer_types ptr() const;
 
    // Convert storage to the logical Cytnx element pointer type.
    //
    // Use ptr_as<T>() for host code and raw memory copies. For complex tensors, T is
    // std::complex<...>. Use gpu_ptr_as<T>() for CUDA kernel code; for complex tensors, T is
    // cuda::std::complex<...>. CUDA library APIs that require cuComplex ABI pointers should cast
    // explicitly at the library-call boundary.
    template <typename T>
    T *ptr_as() const {
      cytnx_error_msg(this->dtype() != Type_class::cy_typeid_v<std::remove_cv_t<T>>,
                      "[ERROR] Attempt to convert dtype %d (%s) to pointer of type %s",
                      this->dtype(), Type_class::getname(this->dtype()).c_str(),
                      Type_class::getname(Type_class::cy_typeid_v<std::remove_cv_t<T>>).c_str());
      return static_cast<T *>(this->_impl->_storage._impl->data());
    }
 
#ifdef UNI_GPU
    // std::variant of pointers to Type_list_gpu, without void ....
    using gpu_pointer_types =
      make_variant_from_transform_t<typename internal::exclude_first<Type_list_gpu>::type,
                                    std::add_pointer>;
 
    // convert this->_impl->_storage->Mem to a typed variant of pointers, excluding void*
    gpu_pointer_types gpu_ptr() const;
 
    // Convert storage to the CUDA kernel element pointer type.
    //
    // This is for CUDA kernel interfaces. Generic host code should use ptr_as<T>() so the logical
    // Cytnx dtype mapping is checked. CUDA library APIs that require cuComplex ABI pointers should
    // cast explicitly at the library-call boundary.
    template <typename T>
    T *gpu_ptr_as() const {
      cytnx_error_msg(
        this->dtype() != Type_class::cy_typeid_gpu_v<std::remove_cv_t<T>>,
        "[ERROR] Attempt to convert dtype %d (%s) to GPU pointer of type %s", this->dtype(),
        Type_class::getname(this->dtype()).c_str(),
        Type_class::getname(Type_class::cy_typeid_gpu_v<std::remove_cv_t<T>>).c_str());
      return static_cast<T *>(this->_impl->_storage._impl->data());
    }
#endif
 
    static Tensor from_storage(const Storage &in) {
      Tensor out;
      boost::intrusive_ptr<Tensor_impl> tmp(new Tensor_impl());
      out._impl = tmp;
      out._impl->Init(in);
      return out;
    }
 
    unsigned int dtype() const { return this->_impl->dtype(); }
 
    int device() const { return this->_impl->device(); }
 
    std::string dtype_str() const { return this->_impl->dtype_str(); }
 
    std::string device_str() const { return this->_impl->device_str(); }
 
    const std::vector<cytnx_uint64> &shape() const { return this->_impl->shape(); }
 
    std::vector<cytnx_int64> strides() const;
 
    cytnx_uint64 rank() const { return this->_impl->rank(); }
 
    cytnx_uint64 size() const { return this->_impl->storage().size(); }
 
    bool is_void() const { return this->_impl->is_void(); }
 
    bool is_scalar() const { return this->_impl->is_scalar(); }
 
    bool is_empty() const { return !this->is_void() && this->size() == 0; }
 
    Tensor clone() const {
      Tensor out;
      out._impl = this->_impl->clone();
      return out;
    }
 
    Tensor to(const int &device) const {
      Tensor out;
      out._impl = this->_impl->to(device);
      return out;
    }
 
    void to_(const int &device) { this->_impl->to_(device); }
 
    const bool &is_contiguous() const { return this->_impl->is_contiguous(); }
 
    Tensor &permute_(const std::vector<cytnx_uint64> &rnks) {
      this->_impl->permute_(rnks);
      return *this;
    }
    template <class... Ts>
    Tensor &permute_(const cytnx_uint64 &e1, const Ts &...elems) {
      std::vector<cytnx_uint64> argv = dynamic_arg_uint64_resolver(e1, elems...);
      this->_impl->permute_(argv);
      return *this;
    }
 
    Tensor permute(const std::vector<cytnx_uint64> &rnks) const {
      Tensor out;
      out._impl = this->_impl->permute(rnks);
      return out;
    }
    template <class... Ts>
    Tensor permute(const cytnx_uint64 &e1, const Ts &...elems) const {
      std::vector<cytnx_uint64> argv = dynamic_arg_uint64_resolver(e1, elems...);
      return this->permute(argv);
    }
 
    Tensor contiguous() const {
      Tensor out;
      out._impl = this->_impl->contiguous();
      return out;
    }
 
    Tensor contiguous_() {
      this->_impl->contiguous_();
      return *this;
    }
 
    Tensor &reshape_(const std::vector<cytnx_int64> &new_shape) {
      this->_impl->reshape_(new_shape);
      return *this;
    }
    Tensor &reshape_(const std::vector<cytnx_uint64> &new_shape) {
      std::vector<cytnx_int64> shape(new_shape.begin(), new_shape.end());
      this->_impl->reshape_(shape);
      return *this;
    }
    Tensor &reshape_(const std::initializer_list<cytnx_int64> &new_shape) {
      std::vector<cytnx_int64> shape = new_shape;
      this->_impl->reshape_(shape);
      return *this;
    }
    template <class... Ts>
    Tensor &reshape_(const cytnx_int64 &e1, const Ts... elems) {
      std::vector<cytnx_int64> shape = dynamic_arg_int64_resolver(e1, elems...);
      this->_impl->reshape_(shape);
      return *this;
    }
 
    Tensor reshape(const std::vector<cytnx_int64> &new_shape) const {
      Tensor out;
      out._impl = this->_impl->reshape(new_shape);
      return out;
    }
 
    Tensor reshape(const std::vector<cytnx_uint64> &new_shape) const {
      std::vector<cytnx_int64> tmp(new_shape.begin(), new_shape.end());
      Tensor out;
      out._impl = this->_impl->reshape(tmp);
      return out;
    }
 
    Tensor reshape(const std::initializer_list<cytnx_int64> &new_shape) const {
      return this->reshape(std::vector<cytnx_int64>(new_shape));
    }
 
    template <class... Ts>
    Tensor reshape(const cytnx_int64 &e1, const Ts &...elems) const {
      std::vector<cytnx_int64> argv = dynamic_arg_int64_resolver(e1, elems...);
      return this->reshape(argv);
    }
 
    Tensor astype(const int &new_type) const {
      Tensor out;
      out._impl = this->_impl->astype(new_type);
      return out;
    }
 
    // Tensor diagonal(){
    //     for(unsigned int i=0;i<this->shape().size();i++){
    //         if(this->shape()[i] != this->shape()[0],"[ERROR] Tensor.diagonal() can only be called
    //         when the subject has equal dimension in each rank.%s","\n");
    //     }
    //
    // }
 
    template <class T>
    T &at(const std::vector<cytnx_uint64> &locator) {
      return this->_impl->at<T>(locator);
    }
 
    template <class T>
    const T &at(const std::vector<cytnx_uint64> &locator) const {
      return this->_impl->at<T>(locator);
    }
    template <class T>
    T &at(const std::initializer_list<cytnx_uint64> &locator) {
      return this->at<T>(std::vector<cytnx_uint64>(locator));
    }
    template <class T>
    const T &at(const std::initializer_list<cytnx_uint64> &locator) const {
      return this->at<T>(std::vector<cytnx_uint64>(locator));
    }
    template <class T, class... Ts>
    const T &at(const cytnx_uint64 &e1, const Ts &...elems) const {
      std::vector<cytnx_uint64> argv = dynamic_arg_uint64_resolver(e1, elems...);
      return this->at<T>(argv);
    }
    template <class T, class... Ts>
    T &at(const cytnx_uint64 &e1, const Ts &...elems) {
      std::vector<cytnx_uint64> argv = dynamic_arg_uint64_resolver(e1, elems...);
      return this->at<T>(argv);
    }
 
    const Scalar::Sproxy at(const std::vector<cytnx_uint64> &locator) const {
      return this->_impl->at(locator);
    }
 
    Scalar::Sproxy at(const std::vector<cytnx_uint64> &locator) { return this->_impl->at(locator); }
    const Scalar::Sproxy at(const std::initializer_list<cytnx_uint64> &locator) const {
      return this->at(std::vector<cytnx_uint64>(locator));
    }
    Scalar::Sproxy at(const std::initializer_list<cytnx_uint64> &locator) {
      return this->at(std::vector<cytnx_uint64>(locator));
    }
 
    template <class T>
    T &item() {
      cytnx_error_msg(this->_impl->storage().size() != 1, "[ERROR][Tensor.item<T>]%s",
                      "item can only be called from a Tensor with only one element\n");
      return this->_impl->storage().at<T>(0);
    }
 
    template <class T>
    const T &item() const {
      cytnx_error_msg(this->_impl->storage().size() != 1, "[ERROR][Tensor.item<T>]%s",
                      "item can only be called from a Tensor with only one element\n");
      return this->_impl->storage().at<T>(0);
    }
 
    const Scalar::Sproxy item() const {
      cytnx_error_msg(this->_impl->storage().size() != 1, "[ERROR][Tensor.item]%s",
                      "item can only be called from a Tensor with only one element\n");
      Scalar::Sproxy out(this->storage()._impl, 0);
      return out;
    }
 
    Scalar::Sproxy item() {
      cytnx_error_msg(this->_impl->storage().size() != 1, "[ERROR][Tensor.item]%s",
                      "item can only be called from a Tensor with only one element\n");
      Scalar::Sproxy out(this->storage()._impl, 0);
      return out;
    }
 
 
    Tensor get(const std::vector<cytnx::Accessor> &accessors,
               std::vector<cytnx_int64> &removed) const {
      Tensor out;
      out._impl = this->_impl->get(accessors, removed);
      return out;
    }
    Tensor get(const std::vector<cytnx::Accessor> &accessors) const {
      Tensor out;
      std::vector<cytnx_int64> removed;
      out._impl = this->_impl->get(accessors, removed);
      return out;
    }
 
    /*
    Tensor get_v2(const std::vector<cytnx::Accessor> &accessors) const{
        Tensor out;
        out._impl = this->_impl->get_v2(accessors);
        return out;
    }
    */
 
    void set(const std::vector<cytnx::Accessor> &accessors, const Tensor &rhs) {
      this->_impl->set(accessors, rhs._impl);
    }
 
    template <class T>
    void set(const std::vector<cytnx::Accessor> &accessors, const T &rc) {
      this->_impl->set(accessors, rc);
    }
    template <class T>
    void set(const std::initializer_list<cytnx::Accessor> &accessors, const T &rc) {
      std::vector<cytnx::Accessor> args = accessors;
      this->set(args, rc);
    }
 
    Storage &storage() const { return this->_impl->storage(); }
 
    template <class T>
    void fill(const T &val) {
      this->_impl->fill(val);
    }
 
    bool equivshape(const Tensor &rhs) {
      if (this->shape() != rhs.shape()) return false;
      return true;
    }
 
    Tensor real();
 
    Tensor imag();
 
    // Arithmic:
    template <class T>
    Tensor &operator+=(const T &rc);
 
    template <class T>
    Tensor &operator-=(const T &rc);
 
    template <class T>
    Tensor &operator*=(const T &rc);
 
    template <class T>
    Tensor &operator/=(const T &rc);
 
    // Tensor &operator+=(const Tproxy &rc);
    // Tensor &operator-=(const Tproxy &rc);
    // Tensor &operator*=(const Tproxy &rc);
    // Tensor &operator/=(const Tproxy &rc);
    /*
    Tensor operator+(const Tproxy &rc){
         return *this + Tensor(rc);
    }
    Tensor operator-(const Tproxy &rc){
         return *this - Tensor(rc);
    }
    Tensor operator*(const Tproxy &rc){
         return *this * Tensor(rc);
    }
    Tensor operator/(const Tproxy &rc){
         return *this / Tensor(rc);
    }
    */
    template <class T>
    Tensor Add(const T &rhs) {
      return *this + rhs;
    }
 
    template <class T>
    Tensor &Add_(const T &rhs) {
      return *this += rhs;
    }
 
    template <class T>
    Tensor Sub(const T &rhs) {
      return *this - rhs;
    }
 
    template <class T>
    Tensor &Sub_(const T &rhs) {
      return *this -= rhs;
    }
 
    template <class T>
    Tensor Mul(const T &rhs) {
      return *this * rhs;
    }
 
    template <class T>
    Tensor &Mul_(const T &rhs) {
      return *this *= rhs;
    }
 
    template <class T>
    Tensor Div(const T &rhs) {
      return *this / rhs;
    }
 
    template <class T>
    Tensor &Div_(const T &rhs) {
      return *this /= rhs;
    }
 
    template <class T>
    Tensor Cpr(const T &rhs) {
      return *this == rhs;
    }
 
    // /**
    //  * @brief Compare each element of the current tensor with the input tensor.
    //  * @details This function Compare each element of the current tensor with the input tensor.
    //  * @param[in] rhs the compared tensor.
    //  */
    // bool approx_eq(const Tensor &rhs, const cytnx_double tol = 0) {
    //   if (this->device() != rhs.device()) {
    //     if (User_debug)
    //       std::cout << "[approx_eq] Tensor device " << this->device()
    //                 << "not equal to rhs tensor device " << rhs.device() << std::endl;
    //     return false;
    //   }
    //   // if (this->dtype() != rhs.dtype()) {
    //   //   std::cout << "[approx_eq] Tensor dtype " << this->dtype()
    //   //             << "not equal to rhs tensor dtype " << rhs.dtype() << std::endl;
    //   //   return false;
    //   // }
    //   if (this->shape() != rhs.shape()) {
    //     if (User_debug)
    //       std::cout << "[approx_eq] Tensor shape " << this->shape()
    //                 << "not equal to rhs tensor shape " << rhs.shape() << std::endl;
    //     return false;
    //   }
    //   if (this->is_contiguous() != rhs.is_contiguous()) {
    //     if (User_debug)
    //       std::cout << "[AreNearlyEqTensor] Tensor contiguous flag " << this->is_contiguous()
    //                 << "not equal to rhs tensor flag " << rhs.is_contiguous() << std::endl;
    //     return false;
    //   }
    //   return this->_impl->_storage.approx_eq(rhs._impl->_storage._impl, tol);
    // }
 
    // template<class T>
    // Tensor& Cpr_(const T &rhs){
    //
    //      return *this == rhs;
    // }
 
    template <class T>
    Tensor Mod(const T &rhs) {
      return *this % rhs;
    }
 
    Tensor operator-() { return this->Mul(-1.); }
 
    Tensor flatten() const {
      Tensor out = this->clone();
      out.contiguous_();
      out.reshape_({-1});
      return out;
    }
 
    void flatten_() {
      this->contiguous_();
      this->reshape_({-1});
    }
 
    void append(const Tensor &rhs) {
      // Tensor in;
      // check Tensor in shape:
      cytnx_error_msg(rhs.is_void() || this->is_void(), "[ERROR] try to append a null Tensor.%s",
                      "\n");
      cytnx_error_msg(this->is_scalar(), "[ERROR] try to append to a rank-0 Tensor.%s", "\n");
      if (!this->is_contiguous()) this->contiguous_();
 
      cytnx_error_msg(rhs.rank() != this->rank() - 1,
                      "[ERROR] try to append a Tensor with rank not match.%s", "\n");
      cytnx_uint64 Nelem = 1;
      for (unsigned int i = 0; i < rhs.shape().size(); i++) {
        cytnx_error_msg(rhs.shape()[i] != this->shape()[i + 1],
                        "[ERROR] dimension mismatch @ rhs.rank: [%d] this: [%d] rhs: [%d]\n", i,
                        this->shape()[i + 1], rhs.shape()[i]);
        Nelem *= rhs.shape()[i];
      }
 
      // check type:
      Tensor in;
      if (rhs.dtype() != this->dtype()) {
        in = rhs.astype(this->dtype());
        if (!in.is_contiguous()) in.contiguous_();
      } else {
        if (!rhs.is_contiguous())
          in = rhs.contiguous();
        else
          in = rhs;
      }
      this->_impl->_shape[0] += 1;
      cytnx_uint64 oldsize = this->_impl->_storage.size();
      this->_impl->_storage.resize(oldsize + Nelem);
      memcpy(((char *)this->_impl->_storage.data()) +
               oldsize * Type.typeSize(this->dtype()) / sizeof(char),
             in._impl->_storage.data(), Type.typeSize(in.dtype()) * Nelem);
    }
    void append(const Storage &srhs) {
      // check Tensor in shape:
      cytnx_error_msg(this->is_void(), "[ERROR] try to append to an uninitialized Tensor.%s", "\n");
      cytnx_error_msg(srhs.size() == 0, "[ERROR] try to append an empty Storage.%s", "\n");
      cytnx_error_msg(this->rank() != 2,
                      "[ERROR] append a storage to Tensor can only accept rank-2 Tensor.%s", "\n");
      if (!this->is_contiguous()) this->contiguous_();
 
      cytnx_error_msg(this->shape().back() != srhs.size(), "[ERROR] Tensor dmension mismatch!%s",
                      "\n");
 
      // check type:
      Storage in;
      if (srhs.dtype() != this->dtype()) {
        in = srhs.astype(this->dtype());
      } else {
        in = srhs;
      }
      this->_impl->_shape[0] += 1;
      cytnx_uint64 oldsize = this->_impl->_storage.size();
      this->_impl->_storage.resize(oldsize + in.size());
      memcpy(((char *)this->_impl->_storage.data()) +
               oldsize * Type.typeSize(this->dtype()) / sizeof(char),
             in._impl->data(), Type.typeSize(in.dtype()) * in.size());
    }
    /*
    void append(const Tensor &rhs){
        // convert to the same type.
        Tensor in;
        if(rhs.dtype() != this->dtype()){
            in = rhs.astype(this->dtype());
        }else{
            in = rhs;
        }
 
        // 1) check rank
        if(this->shape().size()==1){
            // check if rhs is a scalar tensor (only one element)
            cytnx_error_msg(!(rhs.shape().size()==1 && rhs.shape()[0]==1),"[ERROR] trying to append
    a scalar into multidimentional Tensor is not allow.\n Only rank-1 Tensor can accept scalar
    append.%s","\n"); this->_impl->_shape[0]+=1; this->_impl->_storage.append(0);
 
        }else{
            cytnx_error_msg(rhs.shape().size() != this->shape().size()-1,"[ERROR] try to append a
    Tensor with rank not match.%s","\n");
 
        }
        cytnx_error_msg(!this->is_contiguous(),"[ERROR] append require the Tensor to be contiguous.
    suggestion: call contiguous() or contiguous_() first.","\n");
    }
    */
    template <class T>
    void append(const T &rhs) {
      cytnx_error_msg(this->is_void(), "[ERROR] try to append to an uninitialized Tensor.%s", "\n");
      cytnx_error_msg(this->rank() != 1,
                      "[ERROR] trying to append a scalar into multidimentional Tensor is not "
                      "allow.\n Only rank-1 Tensor can accept scalar append.%s",
                      "\n");
      cytnx_error_msg(!this->is_contiguous(),
                      "[ERROR] append require the Tensor to be contiguous. suggestion: call "
                      "contiguous() or contiguous_() first.",
                      "\n");
      this->_impl->_shape[0] += 1;
      this->_impl->_storage.append(rhs);
    }
 
    bool same_data(const Tensor &rhs) const;
 
    // linalg:
    std::vector<Tensor> Svd(const bool &is_UvT = true) const;
 
    std::vector<Tensor> Eigh(const bool &is_V = true, const bool &row_v = false) const;
 
    Tensor &InvM_();
 
    Tensor InvM() const;
 
    Tensor &Inv_(const double &clip = -1.);
 
    Tensor Inv(const double &clip = -1.) const;
 
    Tensor &Conj_();
 
    Tensor Conj() const;
 
    Tensor &Exp_();
 
    Tensor Exp() const;
 
    [[deprecated("use norm() (returns Scalar) instead")]] Tensor Norm() const;
 
    Scalar norm() const;
 
    Tensor Pow(const cytnx_double &p) const;
 
    Tensor &Pow_(const cytnx_double &p);
 
    Tensor Trace(const cytnx_uint64 &a = 0, const cytnx_uint64 &b = 1) const;
 
    Tensor Abs() const;
 
    Tensor &Abs_();
 
    Tensor Max() const;
 
    Tensor Min() const;
 
  };  // class Tensor
 
  // Scalar-like operand concepts for the free Tensor / UniTensor arithmetic and comparison
  // operators. They live here, after the class, because cytnx_scalar_like names Tensor::Tproxy
  // and so needs Tensor to be complete; linalg.hpp includes this header and reuses them.
  //
  // Each concept admits exactly the operand types the corresponding operators are instantiated
  // for, so an unsupported operand is rejected at overload resolution instead of surviving to a
  // link error. Widening one without adding the matching explicit instantiation reintroduces
  // precisely that failure (#1003, operator hygiene -- Ian's review).
 
  // Plain scalar values: the cytnx dtype scalars and cytnx::Scalar. This is the operand set of
  // `operator%` (Tensor and UniTensor) and `operator==` (Tensor), which have no element-proxy
  // instantiations -- see Mod.cpp and Cpr.cpp.
  template <class T>
  concept cytnx_scalar_value =
    CytnxType<std::remove_cvref_t<T>> || std::is_same_v<std::remove_cvref_t<T>, Scalar>;
 
  // Adds both element proxies: the operand set of the free Tensor `+ - * /` operators.
  template <class T>
  concept cytnx_scalar_like =
    cytnx_scalar_value<T> || std::is_same_v<std::remove_cvref_t<T>, Tensor::Tproxy> ||
    std::is_same_v<std::remove_cvref_t<T>, Scalar::Sproxy>;
 
  // Adds only the Scalar proxy: the operand set of the free UniTensor `+ - * /` operators, which
  // are instantiated for Scalar::Sproxy but never for Tensor::Tproxy -- see Add.cpp.
  template <class T>
  concept cytnx_unitensor_scalar_like =
    cytnx_scalar_value<T> || std::is_same_v<std::remove_cvref_t<T>, Scalar::Sproxy>;
 
  // [Note] these are fwd from linalg.hpp; the constraint must match linalg.hpp exactly, or the
  // two declarations become distinct overloads and the unconstrained one stays viable.
  template <cytnx_scalar_like T>
  Tensor operator+(const Tensor &lhs, const T &rc);
  template <cytnx_scalar_like T>
  Tensor operator-(const Tensor &lhs, const T &rhs);
  template <cytnx_scalar_like T>
  Tensor operator*(const Tensor &lhs, const T &rhs);
  template <cytnx_scalar_like T>
  Tensor operator/(const Tensor &lhs, const T &rhs);
 
  Tensor operator+(const Tensor &lhs, const Tensor::Tproxy &rhs);
  Tensor operator-(const Tensor &lhs, const Tensor::Tproxy &rhs);
  Tensor operator*(const Tensor &lhs, const Tensor::Tproxy &rhs);
  Tensor operator/(const Tensor &lhs, const Tensor::Tproxy &rhs);
 
  Tensor operator+(const Tensor &lhs, const Scalar::Sproxy &rhs);
  Tensor operator-(const Tensor &lhs, const Scalar::Sproxy &rhs);
  Tensor operator*(const Tensor &lhs, const Scalar::Sproxy &rhs);
  Tensor operator/(const Tensor &lhs, const Scalar::Sproxy &rhs);
 
  std::ostream &operator<<(std::ostream &os, const Tensor &in);
  std::ostream &operator<<(std::ostream &os, const Tensor::Tproxy &in);
  //{ os << Tensor(in);};
}  // namespace cytnx
 
#endif  // CYTNX_TENSOR_H_