calculate_pTheta.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef cMHN_COMMON_CALCULATE_PTHETA_H
#define cMHN_COMMON_CALCULATE_PTHETA_H
 
#include <cmhn/nontt/jacobi.hpp>
#include <cmhn/nontt/mhn_operator.hpp>
#include <cmhn/tt/als.hpp>
#include <cmhn/tt/mamen.hpp>
#include <cmhn/tt/mhn_operator.hpp>
#include <cmhn/tt/utility.hpp>
 
#include <prc.hpp>
 
namespace cMHN
{
    template<class T, pRC::Size D>
    decltype(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
        [](auto const... Ns)
        {
            return pRC::Tensor<T, Ns...>{};
        }))
    calculatePTheta(nonTT::MHNOperator<T, D> const &op)
    {
        auto const p0 =
            eval(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
                [](auto const... Ns)
                {
                    return pRC::Tensor<T, Ns...>::Single(pRC::identity<T>(),
                        pRC::Subscripts<Ns...>(0));
                }));
 
        return nonTT::jacobi(op, p0);
    }
 
    template<pRC::Size R, class T, class X, pRC::Size D>
    X calculatePTheta(TT::MHNOperator<T, D> const &op, X const &pInit,
        T const &toleranceSolver)
    {
        using ModeSizes = decltype(TT::getModeSizes<D>());
 
        auto const p0 = expand(pRC::makeConstantSequence<pRC::Index, D, 0>(),
            [](auto const... seq)
            {
                return pRC::TensorTrain::Tensor<T, ModeSizes,
                    decltype(pRC::Sizes(pRC::makeConstantSequence<pRC::Size,
                        D - 1, 1>()))>::Single(pRC::identity<T>(), seq...);
            });
 
        auto pTheta = TT::ALS<R>(op, p0, toleranceSolver, eval(pInit));
 
        pTheta /= scalarProduct(pTheta,
            pRC::unit<pRC::TensorTrain::Tensor<T, ModeSizes>>())();
 
        return pTheta;
    }
 
    template<pRC::Size R, class T, pRC::Size D>
    auto calculatePTheta(TT::MHNOperator<T, D> const &op,
        T const &toleranceSolver)
    {
        using ModeSizes = decltype(TT::getModeSizes<D>());
        using Ranks = decltype(TT::getRanks<D, R>());
 
        // use a random pInit
        pRC::SeedSequence seq(8, 16);
        pRC::RandomEngine rng(seq);
        pRC::GaussianDistribution<pRC::Float<>> dist;
        auto pInit = round<Ranks>(
            pRC::random<pRC::TensorTrain::Tensor<T, ModeSizes, Ranks>>(rng,
                dist));
        pInit /= scalarProduct(pInit,
            pRC::unit<pRC::TensorTrain::Tensor<T, ModeSizes>>())();
 
        return calculatePTheta<R>(op, pInit, toleranceSolver);
    }
}
 
#endif // cMHN_COMMON_CALCULATE_PTHETA_H