CoordinateHitAndRunProposal.hpp

Go to the documentation of this file.

#ifndef HOPS_COORDINATEHITANDRUNPROPOSAL_HPP
#define HOPS_COORDINATEHITANDRUNPROPOSAL_HPP
 
#include "ChordStepDistributions.hpp"
#include "../IsSetStepSizeAvailable.hpp"
#include "../../RandomNumberGenerator/RandomNumberGenerator.hpp"
#include <random>
#include "../../FileWriter/CsvWriter.hpp"
 
namespace hops {
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution = UniformStepDistribution<typename MatrixType::Scalar>>
    class CoordinateHitAndRunProposal {
    public:
        using StateType = VectorType;
 
        CoordinateHitAndRunProposal(MatrixType A, VectorType b, VectorType currentState);
 
        void propose(RandomNumberGenerator &randomNumberGenerator);
 
        void acceptProposal();
 
        StateType getState() const;
 
        StateType getProposal() const;
 
        void setState(StateType newState);
 
        void setStepSize(typename MatrixType::Scalar stepSize);
 
        typename MatrixType::Scalar getStepSize() const;
 
        [[nodiscard]] typename MatrixType::Scalar computeLogAcceptanceProbability() {
            return chordStepDistribution.computeInverseNormalizationConstant(0, backwardDistance, forwardDistance)
                   - chordStepDistribution.computeInverseNormalizationConstant(0, backwardDistance - step,
                                                                                 forwardDistance - step);
        }
 
        std::string getName();
 
    private:
        MatrixType A;
        VectorType b;
        StateType state;
        VectorType slacks;
        VectorType inverseDistances;
 
        long coordinateToUpdate = 0;
        typename MatrixType::Scalar step = 0;
        ChordStepDistribution chordStepDistribution;
        typename MatrixType::Scalar forwardDistance;
        typename MatrixType::Scalar backwardDistance;
    };
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::CoordinateHitAndRunProposal(
            MatrixType A_,
            VectorType b_,
            VectorType currentState_) :
            A(std::move(A_)),
            b(std::move(b_)),
            state(std::move(currentState_)) {
        slacks = this->b - this->A * this->state;
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    void CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::propose(
            RandomNumberGenerator &randomNumberGenerator) {
        ++coordinateToUpdate %= state.rows();
 
        inverseDistances = A.col(coordinateToUpdate).cwiseQuotient(slacks);
        forwardDistance = 1. / inverseDistances.maxCoeff();
        backwardDistance = 1. / inverseDistances.minCoeff();
        assert(backwardDistance < 0 && forwardDistance > 0);
        assert(((b - A * state).array() > 0).all());
 
        step = chordStepDistribution.draw(randomNumberGenerator, backwardDistance, forwardDistance);
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    void
    CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::acceptProposal() {
        state(coordinateToUpdate) += step;
        slacks.noalias() -= A.col(coordinateToUpdate) * step;
        step = 0;
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    typename CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::StateType
    CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::getState() const {
        return state;
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    typename CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::StateType
    CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::getProposal() const {
        StateType m_proposal = state;
        m_proposal(coordinateToUpdate) += step;
        return m_proposal;
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    void CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::setState(VectorType newState) {
        CoordinateHitAndRunProposal::state = std::move(newState);
        slacks = b - A * CoordinateHitAndRunProposal::state;
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    void CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::setStepSize(
            typename MatrixType::Scalar stepSize) {
        if constexpr (IsSetStepSizeAvailable<ChordStepDistribution>::value) {
            chordStepDistribution.setStepSize(stepSize);
        } else {
            (void)stepSize;
            throw std::runtime_error("Step size not available.");
        }
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    typename MatrixType::Scalar
    CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::getStepSize() const {
        if constexpr (IsSetStepSizeAvailable<ChordStepDistribution>::value) {
            return chordStepDistribution.getStepSize();
        }
        throw std::runtime_error("Step size not available.");
    }
 
    template<typename MatrixType, typename VectorType, typename ChordStepDistribution>
    std::string CoordinateHitAndRunProposal<MatrixType, VectorType, ChordStepDistribution>::getName() {
        return "Coordinate Hit-and-Run";
    }
}
 
#endif //HOPS_COORDINATEHITANDRUNPROPOSAL_HPP