#pragma once
#include "mps.hpp"
#include <eigen3/Eigen/Eigen>
#include <eigen3/unsupported/Eigen/CXX11/Tensor>
#include <initializer_list>
#include <stdexcept>
namespace e = Eigen;

using namespace std::complex_literals;

/**
 * @brief Construct a matrix product operator in tensor form
 * @matrices N x N initializer lists containing localDim x localDim matrices
 * @details
 *  Constructs a tensor with indices r c i j, where r and c specify the matrix and i j are indices into the matrix
 */
template<typename T>
e::Tensor<T, 4> constructMPO(const std::initializer_list<std::initializer_list<e::MatrixX<T>>>& matrices);

template<typename T>
e::MatrixX<T> MPO_Tensor2Matrix(const e::Tensor<T, 4>& mpo);


/**
 * @brief Reshape a N*localDim x N*localDim matrix into a N x localDim x N x localDim tensor
 */
template<typename T>
e::Tensor<T, 4> bondMatrix2Tensor(const e::MatrixX<T>& bondMatrix, unsigned localDim);


template<typename T>
class Model1D {
    public:
        Model1D(size_t L, unsigned localDim) : L(L), localDim(localDim) {
            if (L <= 2) {  // needs two bonds to function correctly
                throw std::invalid_argument("L must be larger than 2");
            }
        }
        virtual ~Model1D() {};
        inline const std::vector<std::shared_ptr<e::MatrixX<T>>>& getH_BondsMatrices() const { return H_BondsMatrices; }
        inline const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& getH_Bonds() const { return H_Bonds; }
        inline const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& getH_MPO() const { return H_MPO; }
        const size_t L;
        const unsigned localDim;
        T energy(const MPS<T>& psi) const;
        virtual std::string format() const = 0;
    protected:
        virtual void initH_Bonds() = 0;
        virtual void initH_MPO() = 0;
        /// Bonds in tensor form, used for `energy()` function
        std::vector<std::shared_ptr<e::Tensor<T, 4>>> H_Bonds;
        /// Bonds in matrix form, used for calculating time evolution
        std::vector<std::shared_ptr<e::MatrixX<T>>> H_BondsMatrices;
        std::vector<std::shared_ptr<e::Tensor<T, 4>>> H_MPO;
};


/**
 * @brief Class of pauli matrices
 * @details
 *  Intendend for static initialization
 */
template<typename T>
class PauliMatrices {
    public:
        PauliMatrices() :
            sigmaX({{0, 1}, {1, 0}}),
            sigmaY({{0if, -1if}, {1if, 0if}}),
            sigmaZ({{1, 0}, {0, -1}}),
            identity({{1, 0}, {0, 1}}) {}
		inline const e::Matrix2<T>& getSigmaX() const { return sigmaX; }
		inline const e::Matrix2<T>& getSigmaY() const { return sigmaY; }
		inline const e::Matrix2<T>& getSigmaZ() const { return sigmaZ; }
		inline const e::Matrix2<T>& getIdentity() const { return identity; }
    private:
        e::Matrix2<T> sigmaX;
        e::Matrix2<T> sigmaY;
        e::Matrix2<T> sigmaZ;
        e::Matrix2<T> identity;
};