operator.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_TENSOR_TRAIN_OPERATOR_OPERATOR_H
#define pRC_TENSOR_TRAIN_OPERATOR_OPERATOR_H
 
#include <prc/core/basic/limits.hpp>
#include <prc/core/basic/range.hpp>
#include <prc/core/basic/sequence.hpp>
#include <prc/core/complex/type_traits.hpp>
#include <prc/core/functors/mul.hpp>
#include <prc/core/tensor/type_traits.hpp>
#include <prc/core/value/type_traits.hpp>
#include <prc/tensor_train/operator/type_traits.hpp>
#include <prc/tensor_train/operator/views/enumerate.hpp>
#include <prc/tensor_train/operator/views/reference.hpp>
 
namespace pRC::TensorTrain
{
    template<class T, class M, class N, class Ranks, class F>
    Operator(OperatorViews::View<T, M, N, Ranks, F> const &)
        -> Operator<T, M, N, Ranks>;
 
    template<class T, Size... Ms, Size... Ns, Size... Rs>
    class Operator<T, Sizes<Ms...>, Sizes<Ns...>, Sizes<Rs...>,
        If<All<IsSatisfied<(sizeof...(Ms) == sizeof...(Ns))>,
            IsSatisfied<(sizeof...(Ns) - 1 == sizeof...(Rs))>,
            Any<IsValue<T>, IsComplex<T>>>>>
    {
    public:
        using M = pRC::Sizes<Ms...>;
        using N = pRC::Sizes<Ns...>;
        using L = pRC::Sizes<(Ms * Ns)...>;
        using Sizes = pRC::Sizes<Ms..., Ns...>;
 
        using SubscriptsM = pRC::Subscripts<Ms...>;
        using SubscriptsN = pRC::Subscripts<Ns...>;
        using SubscriptsL = pRC::Subscripts<(Ms * Ns)...>;
        using Subscripts = pRC::Subscripts<Ms..., Ns...>;
 
        using Dimension = typename N::Dimension;
 
        using Ranks = pRC::Sizes<Rs...>;
        template<class S,
            If<IsSame<typename S::Dimension, typename Ranks::Dimension>> = 0>
        using ChangeRanks = Operator<T, M, N, S>;
 
        template<Index C>
        using Cores = pRC::Tensor<T, pRC::Sizes<1, Rs..., 1>::size(C),
            M::size(C), N::size(C), pRC::Sizes<1, Rs..., 1>::size(C + 1)>;
 
        using Type = T;
        template<class C>
        using ChangeType = Operator<C, M, N, Ranks>;
 
        using Value = typename T::Value;
        template<class V, If<IsValue<V>> = 0>
        using ChangeValue =
            Operator<typename T::template ChangeValue<V>, M, N, Ranks>;
 
        using Signed = typename T::Signed;
        template<Bool R>
        using ChangeSigned =
            Operator<typename T::template ChangeSigned<R>, M, N, Ranks>;
 
        using Width = typename T::Width;
        template<Size Q>
        using ChangeWidth =
            Operator<typename T::template ChangeWidth<Q>, M, N, Ranks>;
 
        using IsComplexified = typename T::IsComplexified;
        using Complexify = Operator<typename T::Complexify, M, N, Ranks>;
        using NonComplex = Operator<typename T::NonComplex, M, N, Ranks>;
 
        template<class E = typename M::IsLinearizable, If<E> = 0>
        static constexpr auto m()
        {
            return M::size();
        }
 
        static constexpr auto m(Index const dimension)
        {
            return M::size(dimension);
        }
 
        template<class E = typename N::IsLinearizable, If<E> = 0>
        static constexpr auto n()
        {
            return N::size();
        }
 
        static constexpr auto n(Index const dimension)
        {
            return N::size(dimension);
        }
 
        template<class E = typename Sizes::IsLinearizable, If<E> = 0>
        static constexpr auto size()
        {
            return Sizes::size();
        }
 
        static constexpr auto size(Index const dimension)
        {
            return Sizes::size(dimension);
        }
 
        template<class X, class... Is, If<IsConstructible<T, X>> = 0,
            If<IsSatisfied<(sizeof...(Is) == 2 * Dimension())>> = 0>
        static inline constexpr auto Single(X &&value, Is const... indices)
        {
            return Single(forward<X>(value), Subscripts(indices...));
        }
 
        template<class X, If<IsConstructible<T, X>> = 0>
        static inline constexpr auto Single(X &&value,
            Subscripts const &subscripts)
        {
            return expand(makeSeries<Index, Dimension{}>(),
                [&value, &subscripts](auto const... seq)
                {
                    return Single(forward<X>(value),
                        SubscriptsM(subscripts[seq]...),
                        SubscriptsN(subscripts[seq + Dimension()]...));
                });
        }
 
        template<class X, If<IsConstructible<T, X>> = 0>
        static inline constexpr auto Single(X &&value, SubscriptsM const &is,
            SubscriptsN const &js)
        {
            auto const f = [is, js, value = T(forward<X>(value))]<Index C>()
            {
                if constexpr(C == 0)
                {
                    return Cores<C>::Single(value, 0, is[C], js[C], 0);
                }
                else
                {
                    return Cores<C>::Single(identity<T>(), 0, is[C], js[C], 0);
                }
            };
            using F = RemoveConstReference<decltype(f)>;
            return TensorTrain::OperatorViews::Enumerate<T, M, N, Ranks, F>(f);
        }
 
    public:
        ~Operator() = default;
        constexpr Operator(Operator const &) = default;
        constexpr Operator(Operator &&) = default;
        constexpr Operator &operator=(Operator const &) & = default;
        constexpr Operator &operator=(Operator &&) & = default;
        constexpr Operator() = default;
 
        template<class X,
            If<IsAssignable<OperatorViews::Reference<T, M, N, Ranks>, X>> = 0>
        constexpr Operator(X &&other)
        {
            *this = forward<X>(other);
        }
 
        template<class X,
            If<IsAssignable<OperatorViews::Reference<T, M, N, Ranks>, X>> = 0>
        constexpr auto &operator=(X &&rhs) &
        {
            view(*this) = forward<X>(rhs);
            return *this;
        }
 
        template<Index C>
        constexpr decltype(auto) core() &&
        {
            return get<C>(move(mCores));
        }
 
        template<Index C>
        constexpr decltype(auto) core() const &&
        {
            return get<C>(move(mCores));
        }
 
        template<Index C>
        constexpr decltype(auto) core() &
        {
            return get<C>(mCores);
        }
 
        template<Index C>
        constexpr decltype(auto) core() const &
        {
            return get<C>(mCores);
        }
 
        template<class... Is, If<All<IsConvertible<Is, Index>...>> = 0,
            If<IsSatisfied<(sizeof...(Is) == 2 * Dimension())>> = 0>
        constexpr decltype(auto) operator()(Is const... indices) const
        {
            return view(*this)(indices...);
        }
 
        constexpr decltype(auto) operator()(Subscripts const &subscripts) const
        {
            return view(*this)(subscripts);
        }
 
        constexpr decltype(auto) operator()(SubscriptsM const &is,
            SubscriptsN const &js) const
        {
            return view(*this)(is, js);
        }
 
        template<class X, If<IsInvocable<Add, Operator &, X>> = 0>
        constexpr auto &operator+=(X &&rhs) &
        {
            return *this = *this + forward<X>(rhs);
        }
 
        template<class X, If<IsInvocable<Sub, Operator &, X>> = 0>
        constexpr auto &operator-=(X &&rhs) &
        {
            return *this = *this - forward<X>(rhs);
        }
 
        template<class X, If<IsInvocable<Mul, X, Operator &>> = 0>
        constexpr auto &applyOnTheLeft(X &&lhs) &
        {
            view(*this).applyOnTheLeft(forward<X>(lhs));
            return *this;
        }
 
        template<class X, If<IsInvocable<Mul, Operator &, X>> = 0>
        constexpr auto &applyOnTheRight(X &&rhs) &
        {
            view(*this).applyOnTheRight(forward<X>(rhs));
            return *this;
        }
 
        template<class X, If<IsInvocable<Mul, Operator &, X>> = 0>
        constexpr auto &operator*=(X &&rhs) &
        {
            view(*this) *= forward<X>(rhs);
            return *this;
        }
 
        template<class X, If<IsInvocable<Div, Operator &, X>> = 0>
        constexpr auto &operator/=(X &&rhs) &
        {
            return *this = *this / forward<X>(rhs);
        }
 
        template<class E = IsSatisfied<((Ms * ... * 1) * (Ns * ... * 1)) <=
                         NumericLimits<Size>::max() &&
                     typename Sizes::IsLinearizable()>,
            If<E> = 0>
        explicit constexpr operator pRC::Tensor<T, Ms..., Ns...>() const
        {
            return pRC::Tensor(view(*this));
        }
 
    private:
        template<Index... seq>
        static constexpr auto coreTypes(Sequence<Index, seq...>)
        {
            return tuple<Cores<seq>...>{};
        }
 
        using CoreTypes = decltype(coreTypes(makeSeries<Index, Dimension{}>()));
 
    private:
        CoreTypes mCores;
    };
}
 
namespace pRC
{
    template<class T, Size... Ms, Size... Ns, class R>
    Tensor(TensorTrain::Operator<T, Sizes<Ms...>, Sizes<Ns...>, R> const &)
        -> Tensor<T, Ms..., Ns...>;
}
#endif // pRC_TENSOR_TRAIN_OPERATOR_OPERATOR_H