learn_theta.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef cMHN_NONTT_LEARN_THETA_H
#define cMHN_NONTT_LEARN_THETA_H
 
#include <iostream>
#include <map>
#include <string>
#include <tuple>
 
#include <cmhn/common/calculate_pTheta.hpp>
#include <cmhn/nontt/mhn_operator.hpp>
#include <cmhn/utility/read_data.hpp>
#include <cmhn/utility/read_header.hpp>
#include <cmhn/utility/write_theta.hpp>
 
#include <prc.hpp>
 
namespace cMHN::nonTT
{
    template<class T, pRC::Size D, class S>
    std::tuple<pRC::Tensor<T, D, D>,
        std::map<std::string, std::string>,
        std::map<std::string, double>>
    learnTheta(pRC::Tensor<T, D, D> const &theta,
        std::string const &header, std::string const &output,
        std::map<S, T> const &pD, cMHN::Score<T> const &Score,
        cMHN::Regulator<T, D> const &Regulator, T const &toleranceOptimizer,
        T const &toleranceSolverP, T const &toleranceSolverQ)
    {
        auto tempTheta = theta;
 
        auto score = pRC::zero<T>();
 
        pRC::Index at_iter = 0;
 
        pRC::Float<64> startTime = pRC::getTimeInSeconds();
 
        std::map<std::string, double> logInfoNumbers{{"Score", score()},
            {"Iterations", at_iter},
            {"Time", pRC::getTimeInSeconds()() - startTime()},
            {"Lambda", Regulator.lambda()()}};
 
        std::map<std::string, std::string> logInfoNames{
            {"Score Name", Score.name()}, {"Regulator Name", Regulator.name()}};
 
        writeTheta(output, header, tempTheta, logInfoNames, logInfoNumbers);
 
        std::cout << "cMHN learning started (nonTT):" << std::endl;
        std::cout << "\tScore Name:\t" << logInfoNames["Score Name"]
                  << std::endl;
        std::cout << "\tRegulator Name:\t" << logInfoNames["Regulator Name"]
                  << std::endl
                  << std::endl;
        
        // use a random initial pInit
        pRC::SeedSequence seq(8, 16);
        pRC::RandomEngine rng(seq);
        pRC::GaussianDistribution<pRC::Float<>> dist;
        auto pInit = eval(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
            [&](auto const... Ns)
            {
                return pRC::random<pRC::Tensor<T, Ns...>>(rng, dist);
            }));
        pInit = pInit / norm(pInit);
 
        tempTheta = pRC::optimize<pRC::Optimizer::BFGS<>>(
            tempTheta,
            [&output, &at_iter, &score, &pD, &Score, &Regulator, &pInit,
                &toleranceOptimizer, &toleranceSolverP,
                &toleranceSolverQ](auto const &tempTheta, auto &g)
            {
                MHNOperator<T, D> op(tempTheta);
 
                auto const pTheta =
                    ::cMHN::calculatePTheta(op, pInit, toleranceSolverP);
 
                // use last pTheta as next pInit
                pInit = pTheta;
 
                g = pRC::zero();
                score = pRC::zero();
                for(auto const &[k, v] : pD)
                {
                    score += Score.pointwiseScore(v, pTheta(k));
                }
 
                // calculate g = dS/dtheta
                auto rhs =
                    eval(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
                        [](auto const... Ns)
                        {
                            return pRC::zero<pRC::Tensor<T, Ns...>>();
                        }));
 
                for(auto const &[k, v] : pD)
                {
                    rhs(k) = Score.pointwiseDSDP(v, pTheta(k));
                }
 
                auto q = pRC::solve<pRC::Solver::GMRES<>,
                    pRC::Operator::Transform::Transpose>(op, rhs,
                    pRC::zero<decltype(rhs)>(), toleranceSolverQ);
 
                // this follows appendix C of Rudi's Thesis
                for(pRC::Index i = 0; i < D; ++i)
                {
                    auto r = hadamardProduct(q,
                        cMHN::nonTT::applyDerivative(op, pTheta, i));
 
                    pRC::range<pRC::Context::CompileTime, D>(
                        [&](auto const j)
                        {
                            if(i == j)
                            {
                                g(i, j) = pRC::reduce<pRC::Add>(r)();
                            }
                            else
                            {
                                g(i, j) =
                                    pRC::reduce<pRC::Add>(pRC::chip<j>(r, 1))();
                            }
                        });
                }
 
                // add regularization
                // signs b/c pRC::Optimmizer::BFGS minimizes the score function
                g = Regulator.grad(tempTheta) - g;
                return Regulator.score(tempTheta) - score;
            },
            [&output, &header, &score, &at_iter, &startTime, &logInfoNames,
                &logInfoNumbers](auto const &tempTheta)
            {
                logInfoNumbers["Iterations"] = at_iter;
                logInfoNumbers["Score"] = score();
                logInfoNumbers["Time"] =
                    pRC::getTimeInSeconds()() - startTime();
 
                std::cout << "cMHN learning in progress (nonTT):" << std::endl;
                std::cout << std::defaultfloat;
                std::cout << "\tIteration:\t" << logInfoNumbers["Iterations"]
                          << std::endl;
                std::cout << std::scientific;
                std::cout << "\tLambda:\t\t" << logInfoNumbers["Lambda"]
                          << std::endl;
                std::cout << "\tScore:\t\t" << logInfoNumbers["Score"]
                          << std::endl;
                std::cout << "\tTime:\t\t" << logInfoNumbers["Time"]
                          << std::endl;
                std::cout << std::defaultfloat;
 
                writeTheta(output, header, tempTheta, logInfoNames,
                    logInfoNumbers);
                at_iter++;
            },
            toleranceOptimizer);
 
        return std::make_tuple(tempTheta, logInfoNames, logInfoNumbers);
    }
}
 
#endif // cMHN_NONTT_LEARN_THETA_H