tensor.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 <map>
#include <limits>
#include <memory>
#include <random>
#include <functional>

#include "basic.h"
#include "misc/containerSupport.h"
#include "misc/fileIO.h"
#include "misc/random.h"

#include "indexedTensor.h"

namespace xerus {
    //Forwad declarations
    class TensorNetwork;
    
    // NOTE these two functions are used in the template functions of Tensor so they have to be declared before the class..
    class Tensor;
    
    void contract(Tensor& _result, const Tensor& _lhs, const bool _lhsTrans, const Tensor& _rhs, const bool _rhsTrans, const size_t _numModes);
    Tensor contract(const Tensor& _lhs, const bool _lhsTrans, const Tensor& _rhs, const bool _rhsTrans, const size_t _numModes);
    
    
    
    void reshuffle(Tensor& _out, const Tensor& _base, const std::vector<size_t>& _shuffle);
    Tensor reshuffle(const Tensor& _base, const std::vector<size_t>& _shuffle);
    
    
    class Tensor final {
    public:
        static size_t sparsityFactor; // NOTE not const so that users can modify this value!
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Auxiliary types- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        enum class Initialisation : bool { None, Zero };
        
        enum class Representation : bool { Dense, Sparse };
        
        typedef std::vector<size_t> DimensionTuple; // NOTE must not be declared as "using.." (internal segfault in gcc4.8.1)
        
        typedef std::vector<size_t> MultiIndex; // NOTE must not be declared as "using.." (internal segfault in gcc4.8.1)
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Member variables - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        DimensionTuple dimensions;
        
        size_t size = 1;
        
        value_t factor = 1.0;
        
        Representation representation = Representation::Sparse;
        
    private:
        std::shared_ptr<value_t> denseData;
        
        std::shared_ptr<std::map<size_t, value_t>> sparseData;
        
    public:
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Constructors - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        explicit Tensor(const Representation _representation = Representation::Sparse);
        
        
        Tensor( const Tensor&  _other ) = default;
        
        
        Tensor(       Tensor&& _other ) noexcept = default;
        
        explicit Tensor(DimensionTuple _dimensions, const Representation _representation = Representation::Sparse, const Initialisation _init = Initialisation::Zero);
        
        
        template<XERUS_ADD_MOVE(Vec, DimensionTuple), XERUS_ADD_MOVE(SPtr, std::shared_ptr<value_t>)>
        explicit Tensor(Vec&& _dimensions, SPtr&& _data)
        : dimensions(std::forward<Vec>(_dimensions)), size(misc::product(dimensions)), representation(Representation::Dense), denseData(std::forward<SPtr>(_data)) { }
        
        
        explicit Tensor(DimensionTuple _dimensions, std::unique_ptr<value_t[]>&& _data);
        
        
        explicit Tensor(DimensionTuple _dimensions, const std::function<value_t()>& _f);
        
        
        explicit Tensor(DimensionTuple _dimensions, const std::function<value_t(const size_t)>& _f);
        
        
        explicit Tensor(DimensionTuple _dimensions, const std::function<value_t(const MultiIndex&)>& _f);
        
        
        Tensor(DimensionTuple _dimensions, const size_t _N, const std::function<std::pair<size_t, value_t>(size_t, size_t)>& _f);
        
        
        template<class distribution=std::normal_distribution<value_t>, class generator=std::mt19937_64>
        static Tensor XERUS_warn_unused random(DimensionTuple _dimensions, distribution& _dist=xerus::misc::defaultNormalDistribution, generator& _rnd=xerus::misc::randomEngine) {
            Tensor result(std::move(_dimensions), Representation::Dense, Initialisation::None);
            value_t* const dataPtr = result.denseData.get();
            for(size_t i = 0; i < result.size; ++i) {
                dataPtr[i] = _dist(_rnd);
            }
            return result;
        }
        
        
        template<class distribution=std::normal_distribution<value_t>, class generator=std::mt19937_64>
        XERUS_force_inline static Tensor XERUS_warn_unused random(std::initializer_list<size_t>&& _dimensions, distribution& _dist=xerus::misc::defaultNormalDistribution, generator& _rnd=xerus::misc::randomEngine) {
            return Tensor::random(DimensionTuple(std::move(_dimensions)), _dist, _rnd);
        }
        
        
        template<class generator=std::mt19937_64>
        static Tensor XERUS_warn_unused random_orthogonal(DimensionTuple _dimensions1, DimensionTuple _dimensions2, generator& _rnd=xerus::misc::randomEngine) {
            std::vector<size_t> dimensions = _dimensions1;
            dimensions.insert(dimensions.end(), _dimensions2.begin(), _dimensions2.end());
            const size_t m = misc::product(_dimensions1);
            const size_t n = misc::product(_dimensions2);
            const size_t max = std::max(m,n);
            const size_t min = std::min(m,n);
            Tensor result({max,min}, Representation::Sparse, Initialisation::Zero);
            
            typename generator::result_type randomness = 0;
            const size_t restart = size_t(std::log2(generator::max()));
            for (size_t i=0; i<min; ++i) {
                auto idx = i%restart;
                if (idx == 0) {
                    randomness = _rnd();
                }
                if (randomness & (1<<idx)) {
                    result[{i,i}] = 1;
                } else {
                    result[{i,i}] = -1;
                }
            }
            
            for (size_t i=0; i<min-1; ++i) { // do k = n-1 to 1 by -1; 
                Tensor u = Tensor::random({max-i}, misc::defaultNormalDistribution, _rnd);
                u[0] -= u.frob_norm();
                u /= u.frob_norm();
                u.apply_factor();
                contract(u, u, false, u, false, 0);
                u *= -2.0;
                Tensor p = Tensor::identity({max,max});
                p.offset_add(u, {i,i});
                contract(result, p, false, result, false, 1);
            }
            
            if (m != max) {
                reshuffle(result, result, {1,0});
            }
            result.reinterpret_dimensions(std::move(dimensions));
            
            return result;
        }
        
        
        template<class distribution=std::normal_distribution<value_t>, class generator=std::mt19937_64>
        static Tensor XERUS_warn_unused random(DimensionTuple _dimensions, const size_t _N, distribution& _dist=xerus::misc::defaultNormalDistribution, generator& _rnd=xerus::misc::randomEngine) {
            Tensor result(std::move(_dimensions), Representation::Sparse, Initialisation::Zero);
            XERUS_REQUIRE(_N <= result.size, " Cannot create " << _N << " non zero entries in a tensor with only " << result.size << " total entries!");
            
            std::uniform_int_distribution<size_t> entryDist(0, result.size-1);
            while(result.sparseData->size() < _N) {
                result.sparseData->emplace(entryDist(_rnd), _dist(_rnd));
            }
            return result;
        }
        
        
        template<class distribution=std::normal_distribution<value_t>, class generator=std::mt19937_64>
        XERUS_force_inline static Tensor XERUS_warn_unused random(std::initializer_list<size_t>&& _dimensions, const size_t _N, distribution& _dist, generator& _rnd) {
            return Tensor::random(DimensionTuple(_dimensions), _N, _rnd, _dist);
        }
        
        
        static Tensor XERUS_warn_unused ones(DimensionTuple _dimensions);
        
        
        static Tensor XERUS_warn_unused identity(DimensionTuple _dimensions);
        
        
        static Tensor XERUS_warn_unused kronecker(DimensionTuple _dimensions);
        
        
        static Tensor XERUS_warn_unused dirac(DimensionTuple _dimensions, const MultiIndex& _position);
        
        
        static Tensor XERUS_warn_unused dirac(DimensionTuple _dimensions, const size_t _position);
        

        Tensor dense_copy() const;
        
        
        Tensor sparse_copy() const;
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Standard operators - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        Tensor& operator=(const Tensor&  _other) = default;
        
        
        Tensor& operator=(      Tensor&& _other) = default;
        
        
        Tensor& operator=(const TensorNetwork& _network);
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Information - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        size_t degree() const;
        
        bool has_factor() const;
        
        bool is_dense() const;
        
        bool is_sparse() const;
        
        size_t sparsity() const;
        
        size_t count_non_zero_entries(const value_t _eps = std::numeric_limits<value_t>::epsilon()) const;
        
        bool all_entries_valid() const;
        
        size_t reorder_cost() const;
        
        value_t frob_norm() const;
        
        value_t one_norm() const;
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Basic arithmetics - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        Tensor& operator+=(const Tensor& _other);
        
        Tensor& operator-=(const Tensor& _other);
        
        Tensor& operator*=(const value_t _factor);
        
        Tensor& operator/=(const value_t _divisor);
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Access - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        value_t& operator[](const size_t _position);
        
        value_t operator[](const size_t _position) const;
        
        value_t& operator[](const MultiIndex& _positions);
        
        value_t operator[](const MultiIndex& _positions) const;
        
        
        value_t& at(const size_t _position);
        
        value_t cat(const size_t _position) const;
        
        value_t* get_dense_data();
        
        value_t* get_unsanitized_dense_data();
        
        const value_t* get_unsanitized_dense_data() const;
        
        value_t* override_dense_data();
        
        const std::shared_ptr<value_t>& get_internal_dense_data();
        
        std::map<size_t, value_t>& get_sparse_data();
        
        std::map<size_t, value_t>& get_unsanitized_sparse_data();
        
        const std::map<size_t, value_t>& get_unsanitized_sparse_data() const;
        
        std::map<size_t, value_t>& override_sparse_data();
        
        const std::shared_ptr<std::map<size_t, value_t>>& get_internal_sparse_data();
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Indexing - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        template<typename... args>
        internal::IndexedTensor<Tensor> operator()(args... _args) {
            return internal::IndexedTensor<Tensor>(this, std::vector<Index>({_args...}), false);
        }
        
        
        template<typename... args>
        internal::IndexedTensorReadOnly<Tensor> operator()(args... _args) const {
            return internal::IndexedTensorReadOnly<Tensor>(this, std::vector<Index>({_args...}));
        }
        
        
        internal::IndexedTensor<Tensor> operator()(const std::vector<Index>&  _indices);
        
        
        internal::IndexedTensor<Tensor> operator()(   std::vector<Index>&& _indices);

        
        internal::IndexedTensorReadOnly<Tensor> operator()(const std::vector<Index>&  _indices) const;
        
        
        internal::IndexedTensorReadOnly<Tensor> operator()(   std::vector<Index>&& _indices) const;
        
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Modifiers - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        void reset(DimensionTuple _newDim, const Representation _representation, const Initialisation _init = Initialisation::Zero);
        
        
        void reset(DimensionTuple _newDim, const Initialisation _init = Initialisation::Zero);
        
        void reset();
        
        
        void reset(DimensionTuple _newDim, const std::shared_ptr<value_t>& _newData);
        
        
        void reset(DimensionTuple _newDim, std::unique_ptr<value_t[]>&& _newData);
        
        
        void reset(DimensionTuple _newDim, std::map<size_t, value_t>&& _newData);
        
        
        void reinterpret_dimensions(DimensionTuple _newDimensions);
        
        void resize_mode(const size_t _mode, const size_t _newDim, size_t _cutPos=~0ul);
        
        
        void fix_mode(const size_t _mode, const size_t _slatePosition);
        
        
        void remove_slate(const size_t _mode, const size_t _pos);
        
        
        void perform_trace(size_t _firstMode, size_t _secondMode);
        
        
        void modify_diagonal_entries(const std::function<void(value_t&)>& _f);
        
        
        void modify_diagonal_entries(const std::function<void(value_t&, const size_t)>& _f);
        
        
        void modify_entries(const std::function<void(value_t&)>& _f);
        
        
        void modify_entries(const std::function<void(value_t&, const size_t)>& _f);
        
        
        void modify_entries(const std::function<void(value_t&, const MultiIndex&)>& _f);
        
        void offset_add(const Tensor& _other, const std::vector<size_t>& _offsets);
        
        void use_dense_representation();
        
        
        void use_dense_representation_if_desirable();
        
        
        void use_sparse_representation(const value_t _eps = std::numeric_limits<value_t>::epsilon());
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Miscellaneous - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        std::string to_string() const;
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Auxiliary functions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
        
        static size_t multiIndex_to_position(const MultiIndex& _multiIndex, const DimensionTuple& _dimensions);
        
        static MultiIndex position_to_multiIndex(size_t _position, const DimensionTuple& _dimensions);
        
        
        /*- - - - - - - - - - - - - - - - - - - - - - - - - - Internal Helper functions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
    protected:
        template<int sign>
        static void plus_minus_equal(Tensor& _me, const Tensor& _other);
        
        static void add_sparse_to_full(const std::shared_ptr<value_t>& _denseData, const value_t _factor, const std::shared_ptr<const std::map<size_t, value_t>>& _sparseData);
        
        static void add_sparse_to_sparse(const std::shared_ptr<std::map<size_t, value_t>>& _sum, const value_t _factor, const std::shared_ptr<const std::map<size_t, value_t>>& _summand);
        
    public:
        
        void ensure_own_data();
        
        void ensure_own_data_no_copy();
        
        void apply_factor();
        
        void ensure_own_data_and_apply_factor();
    };
    
    /*- - - - - - - - - - - - - - - - - - - - - - - - - - External functions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
    
    XERUS_force_inline void contract(Tensor& _result, const Tensor& _lhs,  const Tensor& _rhs, const size_t _numModes) {
        contract(_result, _lhs, false, _rhs, false, _numModes);
    }
    
    XERUS_force_inline Tensor contract(const Tensor& _lhs, const Tensor& _rhs, const size_t _numModes) {
        return contract(_lhs, false, _rhs, false, _numModes);
    }
    
    /*- - - - - - - - - - - - - - - - - - - - - - - - - - Basic arithmetics - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
    
    Tensor operator+(Tensor _lhs, const Tensor& _rhs);
    
    Tensor operator-(Tensor _lhs, const Tensor& _rhs);
    
    Tensor operator*(const value_t _factor, Tensor _tensor);
    
    Tensor operator*(Tensor _tensor, const value_t _factor);
    
    Tensor operator/(Tensor _tensor, const value_t _divisor);
    
    
    
    static XERUS_force_inline value_t frob_norm(const Tensor& _tensor) { return _tensor.frob_norm(); }
    
    static XERUS_force_inline value_t one_norm(const Tensor& _tensor) { return _tensor.one_norm(); }
    
    void calculate_svd(Tensor& _U, Tensor& _S, Tensor& _Vt, Tensor _input, const size_t _splitPos, const size_t _maxRank, const value_t _eps);
    
    void calculate_qr(Tensor& _Q, Tensor& _R, Tensor _input, const size_t _splitPos);
    
    void calculate_rq(Tensor& _R, Tensor& _Q, Tensor _input, const size_t _splitPos);
    
    void calculate_qc(Tensor& _Q, Tensor& _C, Tensor _input, const size_t _splitPos);
    
    void calculate_cq(Tensor& _C, Tensor& _Q, Tensor _input, const size_t _splitPos);
    
    void pseudo_inverse(Tensor& _inverse, const Tensor& _input, const size_t _splitPos);
    Tensor pseudo_inverse(const Tensor& _input, const size_t _splitPos);
    
    void solve_least_squares(Tensor& _X, const Tensor& _A, const Tensor& _B, const size_t _extraDegree = 0);
    
    void solve(Tensor &_X, const Tensor &_A, const Tensor &_B, size_t _extraDegree = 0);
    
    Tensor entrywise_product(const Tensor &_A, const Tensor &_B);
    
    bool approx_equal(const xerus::Tensor& _a, const xerus::Tensor& _b, const xerus::value_t _eps = EPSILON);
    
    bool approx_entrywise_equal(const xerus::Tensor& _a, const xerus::Tensor& _b, const xerus::value_t _eps = EPSILON);
    
    
    bool approx_entrywise_equal(const xerus::Tensor& _tensor, const std::vector<value_t>& _values, const xerus::value_t _eps = EPSILON);
    
    std::ostream& operator<<(std::ostream& _out, const Tensor& _tensor);
    
    namespace misc {
        void stream_writer(std::ostream &_stream, const Tensor &_obj, const FileFormat _format);

        void stream_reader(std::istream &_stream, Tensor &_obj, const FileFormat _format);
    }
}