ttNetwork.h

Go to the documentation of this file.

// Xerus - A General Purpose Tensor Library
// Copyright (C) 2014-2017 Benjamin Huber and Sebastian Wolf. 
// 
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
// 
// Xerus 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 Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org 
// or contact us at contact@libXerus.org.

#pragma once

 
#include "misc/containerSupport.h"
#include "misc/check.h"

#include "index.h"
#include "indexedTensor_tensor_factorisations.h"
#include "tensor.h"
#include "tensorNetwork.h"
#include "indexedTensor.h"
#include "indexedTensorMoveable.h"
#include "indexedTensorList.h"

namespace xerus {
    template<bool isOperator>
    class TTNetwork final : public TensorNetwork {
    public:
        static constexpr const size_t N = isOperator?2:1;
        
        bool canonicalized;
        
        size_t corePosition;        
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Constructors - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        explicit TTNetwork();
        
        
        TTNetwork(const TTNetwork & _cpy) = default;
        
        
        TTNetwork(      TTNetwork&& _mov) noexcept = default;
        
        
        explicit TTNetwork(const size_t _degree);
        
        
        explicit TTNetwork(Tensor::DimensionTuple _dimensions);
        
        
        explicit TTNetwork(const Tensor& _tensor, const double _eps=EPSILON, const size_t _maxRank=std::numeric_limits<size_t>::max());
        
        
        explicit TTNetwork(const Tensor& _tensor, const double _eps, const RankTuple& _maxRanks);
        
        
//      explicit TTNetwork(const TensorNetwork &_network, double _eps=EPSILON);
        
        
        template<class distribution=std::normal_distribution<value_t>, class generator=std::mt19937_64>
        static TTNetwork XERUS_warn_unused random(std::vector<size_t> _dimensions, const std::vector<size_t> &_ranks, distribution& _dist=xerus::misc::defaultNormalDistribution, generator& _rnd=xerus::misc::randomEngine) {
            const size_t numComponents = _dimensions.size()/N;
            XERUS_REQUIRE(_dimensions.size()%N==0, "Illegal number of dimensions for TTOperator.");
            XERUS_REQUIRE(_ranks.size()+1 == numComponents,"Non-matching amount of ranks given to TTNetwork::random.");
            XERUS_REQUIRE(!misc::contains(_dimensions, size_t(0)), "Trying to construct a TTTensor with dimension 0 is not possible.");
            XERUS_REQUIRE(!misc::contains(_ranks, size_t(0)), "Trying to construct random TTTensor with rank 0 is illegal.");
            
            
            
            TTNetwork result(_dimensions.size());
            const std::vector<size_t> targetRank = reduce_to_maximal_ranks(_ranks, _dimensions);
            
            for(size_t i = 0; i < numComponents; ++i) {
                const size_t leftRank = i==0 ? 1 : targetRank[i-1];
                const size_t rightRank = (i==numComponents-1) ? 1 : targetRank[i];

                if(isOperator) {
                    const auto rndComp = Tensor::random({leftRank, _dimensions[i], _dimensions[numComponents+i], rightRank}, _dist, _rnd);
                    result.set_component(i, rndComp);
                } else {
                    const auto rndComp = Tensor::random({leftRank, _dimensions[i], rightRank}, _dist, _rnd);
                    result.set_component(i, rndComp);
                }
            }
            result.move_core(0);
            return result;
        }
        
        
        template<class distribution=std::normal_distribution<value_t>, class generator=std::mt19937_64>
        static TTNetwork XERUS_warn_unused random(const std::vector<size_t>& _dimensions, const size_t _rank, distribution& _dist=xerus::misc::defaultNormalDistribution, generator& _rnd=xerus::misc::randomEngine) {
            return TTNetwork::random(_dimensions, std::vector<size_t>(_dimensions.size()/N-1, _rank), _dist, _rnd);
        }
        
        
        template<class distribution=std::normal_distribution<value_t>, class generator=std::mt19937_64>
        static TTNetwork XERUS_warn_unused random(const std::vector<size_t>& _dimensions, const std::vector<size_t> &_ranks, 
                                const std::function<void(Tensor&)> &_modifySingularValues,
                                distribution& _dist=xerus::misc::defaultNormalDistribution, generator& _rnd=xerus::misc::randomEngine) 
        {
            const size_t numComponents = _dimensions.size()/N;
            
            TTNetwork result = random(_dimensions, _ranks, _dist, _rnd);
            
            const Index i,j,k,l,m;
            
            for (size_t pos = 0; pos+1 < numComponents; ++pos) {
                Tensor A;
                A(i,j^N,k^N,l) = result.component(pos)(i,j^N,m) * result.component(pos+1)(m,k^N,l);
                Tensor U,S,Vt;
                (U(i&1,j),S(j,k),Vt(k,l&1)) = SVD(A(i/2,l/2));
                _modifySingularValues(S);
                Vt(j,l&1) = S(j,k) * Vt(k,l&1);
                result.set_component(pos, U);
                result.set_component(pos+1, Vt);
                result.assume_core_position(pos+1);
            }
            
            result.require_correct_format();
            XERUS_INTERNAL_CHECK(!result.exceeds_maximal_ranks(), "Internal Error");
            result.canonicalize_left();
            return result;
        }
        
        
        static TTNetwork XERUS_warn_unused ones(const std::vector<size_t>& _dimensions);
        
        
        template<bool B = isOperator, typename std::enable_if<B, int>::type = 0>
        static TTNetwork XERUS_warn_unused identity(const std::vector<size_t>& _dimensions);
        
        
        static TTNetwork XERUS_warn_unused kronecker(const std::vector<size_t>& _dimensions);
        
        
        static TTNetwork XERUS_warn_unused dirac(std::vector<size_t> _dimensions, const std::vector<size_t>& _position);
        
        
        static TTNetwork XERUS_warn_unused dirac(std::vector<size_t> _dimensions, const size_t _position);
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Standard Operators - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        TTNetwork& operator=(const TTNetwork&  _other) = default;
        
        
        TTNetwork& operator=(      TTNetwork&& _other) = default;
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Internal helper functions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
    protected:
        static void construct_train_from_full(TensorNetwork& _out, const Tensor& _A, const double _eps);
        
        
        bool exceeds_maximal_ranks() const;
        
        
        size_t num_ranks() const;
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Miscellaneous - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
    public:
        static std::vector<size_t> reduce_to_maximal_ranks(std::vector<size_t> _ranks, const std::vector<size_t>& _dimensions);
        
        static size_t degrees_of_freedom(const std::vector<size_t> &_dimensions, const std::vector<size_t> &_ranks);
        
        size_t degrees_of_freedom();
        
        virtual void fix_mode(const size_t _mode, const size_t _slatePosition) override;
        
        virtual void resize_mode(const size_t _mode, const size_t _newDim, const size_t _cutPos=~0ul) override;
        
        void use_dense_representations();
        
        Tensor& component(const size_t _idx);
        
        
        const Tensor& get_component(const size_t _idx) const;
        
        
        void set_component(const size_t _idx, Tensor _T);
        
        
        std::pair<TensorNetwork, TensorNetwork> chop(const size_t _position) const;
        
        
        void round(const std::vector<size_t>& _maxRanks, const double _eps = EPSILON);
        
        
        void round(const size_t _maxRank);
        
        
        void round(const int _maxRank);
        
        
        void round(const value_t _eps);
        
        
        void soft_threshold(const double _tau, const bool _preventZero = false);
        
        
        void soft_threshold(const std::vector<double>& _taus, const bool _preventZero = false);
        
        
        std::vector<size_t> ranks() const;
        
        
        size_t rank(const size_t _i) const;
        
        
        void move_core(const size_t _position, const bool _keepRank=false);
        
        
        void assume_core_position(const size_t _pos);
        
        
        void canonicalize_left();
        
        
        void canonicalize_right();
        
        
        template<bool B = isOperator, typename std::enable_if<B, int>::type = 0>
        void transpose() {
            const std::vector<size_t> shuffle({0,2,1,3});
            for (size_t n = 0; n < degree()/N; ++n) {
                xerus::reshuffle(component(n), component(n), shuffle);
            }
        }
        
        
        virtual TensorNetwork* get_copy() const override;
        
        
        virtual void contract_unconnected_subnetworks() override;
        
        
        virtual value_t frob_norm() const override;
        
        
        virtual void require_correct_format() const override;
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - -  Basic arithmetics - - - - - - - - - - - - - - - - - - - - - - - - - - */
        TTNetwork& operator+=(const TTNetwork& _other);
        
        
        TTNetwork& operator-=(const TTNetwork& _other);
        
        virtual void operator*=(const value_t _factor) override;
        
        virtual void operator/=(const value_t _divisor) override;
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Operator specializations - - - - - - - - - - - - - - - - - - - - - - - - - - */
        static bool specialized_contraction_f(std::unique_ptr<internal::IndexedTensorMoveable<TensorNetwork>>& _out, internal::IndexedTensorReadOnly<TensorNetwork>&& _me, internal::IndexedTensorReadOnly<TensorNetwork>&& _other);
        
        static bool specialized_sum_f(std::unique_ptr<internal::IndexedTensorMoveable<TensorNetwork>>& _out, internal::IndexedTensorReadOnly<TensorNetwork>&& _me, internal::IndexedTensorReadOnly<TensorNetwork>&& _other);
        
        virtual bool specialized_contraction(std::unique_ptr<internal::IndexedTensorMoveable<TensorNetwork>>& _out, internal::IndexedTensorReadOnly<TensorNetwork>&& _me, internal::IndexedTensorReadOnly<TensorNetwork>&& _other) const override {
            return specialized_contraction_f(_out, std::move(_me), std::move(_other));
        }
        
        virtual bool specialized_sum(std::unique_ptr<internal::IndexedTensorMoveable<TensorNetwork>>& _out, internal::IndexedTensorReadOnly<TensorNetwork>&& _me, internal::IndexedTensorReadOnly<TensorNetwork>&& _other) const override {
            return specialized_sum_f(_out, std::move(_me), std::move(_other));
        }
        
        virtual void specialized_evaluation(internal::IndexedTensorWritable<TensorNetwork>&& _me, internal::IndexedTensorReadOnly<TensorNetwork>&& _other) override;
        
    };

    typedef TTNetwork<false> TTTensor;
    typedef TTNetwork<true> TTOperator;
    

    template<bool isOperator>
    TTNetwork<isOperator> operator+(TTNetwork<isOperator> _lhs, const TTNetwork<isOperator>& _rhs);


    template<bool isOperator>
    TTNetwork<isOperator> operator-(TTNetwork<isOperator> _lhs, const TTNetwork<isOperator>& _rhs);
    
    
    template<bool isOperator>
    TTNetwork<isOperator> operator*(TTNetwork<isOperator> _network, const value_t _factor);
    
    
    template<bool isOperator>
    TTNetwork<isOperator> operator*(const value_t _factor, TTNetwork<isOperator> _network);


    template<bool isOperator>
    TTNetwork<isOperator> operator/(TTNetwork<isOperator> _network, const value_t _divisor);
    
    
    template<bool isOperator>
    TTNetwork<isOperator> entrywise_product(const TTNetwork<isOperator>& _A, const TTNetwork<isOperator>& _B);
    
    template<bool isOperator>
    TTNetwork<isOperator> dyadic_product(const TTNetwork<isOperator> &_lhs, const TTNetwork<isOperator> &_rhs);
    
    
    template<bool isOperator>
    TTNetwork<isOperator> dyadic_product(const std::vector<TTNetwork<isOperator>> &_tensors);
    
    
    namespace misc {
        
        template<bool isOperator>
        void stream_writer(std::ostream& _stream, const TTNetwork<isOperator> &_obj, misc::FileFormat _format);
        
        template<bool isOperator>
        void stream_reader(std::istream& _stream, TTNetwork<isOperator> &_obj, const misc::FileFormat _format);
    }
}