Value Iteration Implementation for MDPs

Question

I've been working for a while on a decision theory library, and since I've never really had any formal training in code best practices I'd love to hear your feedback. This particular class is one of my older ones, and it performs the Value Iteration algorithm on a supplied class which respects a particular interface.

There are includes in the code, but I'm not sure I should include all relevant code, or give more information about parts not shown. Please tell me if I should.

My main concern is whether I should live getters and setters like this, or if I should templatize/make every method inline so as to improve ease of use/speed. Of course, if there is anything else you feel I should change, I'd be very happy to know.

Header file:

#ifndef AI_TOOLBOX_MDP_VALUE_ITERATION_HEADER_FILE
#define AI_TOOLBOX_MDP_VALUE_ITERATION_HEADER_FILE

#include <tuple>
#include <iostream>
#include <iterator>

#include <AIToolbox/MDP/Types.hpp>
#include <AIToolbox/MDP/Utils.hpp>
#include <AIToolbox/ProbabilityUtils.hpp>

namespace AIToolbox {
    namespace MDP {
        /**
         * @brief This class applies the value iteration algorithm on a Model.
         *
         * This algorithm solves an MDP model for the specified horizon, or less
         * if convergence is encountered.
         *
         * The idea of this algorithm is to iteratively compute the
         * ValueFunction for the MDP optimal policy. On the first iteration,
         * the ValueFunction for horizon 1 is obtained. On the second
         * iteration, the one for horizon 2. This process is repeated until the
         * ValueFunction has converged to a specific value within a certain
         * accuracy, or the horizon requested is reached.
         *
         * This implementation in particular is ported from the MATLAB
         * MDPToolbox (although it is simplified).
         */
        class ValueIteration {
            public:
                /**
                 * @brief Basic constructor.
                 *
                 * The epsilon parameter must be >= 0.0, otherwise the
                 * constructor will throw an std::runtime_error. The epsilon
                 * parameter sets the convergence criterion. An epsilon of 0.0
                 * forces ValueIteration to perform a number of iterations
                 * equal to the horizon specified. Otherwise, ValueIteration
                 * will stop as soon as the difference between two iterations
                 * is less than the epsilon specified.
                 *
                 * Note that the default value function size needs to match
                 * the number of states of the Model. Otherwise it will
                 * be ignored. An empty value function will be defaulted
                 * to all zeroes.
                 *
                 * @param horizon The maximum number of iterations to perform.
                 * @param epsilon The epsilon factor to stop the value iteration loop.
                 * @param v The initial value function from which to start the loop.
                 */
                ValueIteration(unsigned horizon, double epsilon = 0.001, ValueFunction v = ValueFunction(Values(0), Actions(0)));

                /**
                 * @brief This function applies value iteration on an MDP to solve it.
                 *
                 * The algorithm is constrained by the currently set parameters.
                 *
                 * @tparam M The type of the solvable MDP.
                 * @param m The MDP that needs to be solved.
                 * @return A tuple containing a boolean value specifying whether
                 *         the specified epsilon bound was reached and the
                 *         ValueFunction and the QFunction for the Model.
                 */
                template <typename M, typename std::enable_if<is_model<M>::value, int>::type = 0>
                std::tuple<bool, ValueFunction, QFunction> operator()(const M & m);

                /**
                 * @brief This function sets the epsilon parameter.
                 *
                 * The epsilon parameter must be >= 0.0, otherwise the
                 * constructor will throw an std::runtime_error. The epsilon
                 * parameter sets the convergence criterion. An epsilon of 0.0
                 * forces ValueIteration to perform a number of iterations
                 * equal to the horizon specified. Otherwise, ValueIteration
                 * will stop as soon as the difference between two iterations
                 * is less than the epsilon specified.
                 *
                 * @param e The new epsilon parameter.
                 */
                void setEpsilon(double e);

                /**
                 * @brief This function sets the horizon parameter.
                 *
                 * @param h The new horizon parameter.
                 */
                void setHorizon(unsigned h);

                /**
                 * @brief This function sets the starting value function.
                 *
                 * An empty value function defaults to all zeroes. Note
                 * that the default value function size needs to match
                 * the number of states of the Model that needs to be
                 * solved. Otherwise it will be ignored.
                 *
                 * @param v The new starting value function.
                 */
                void setValueFunction(ValueFunction v);

                /**
                 * @brief This function will return the currently set epsilon parameter.
                 *
                 * @return The currently set epsilon parameter.
                 */
                double getEpsilon() const;

                /**
                 * @brief This function will return the current horizon parameter.
                 *
                 * @return The currently set horizon parameter.
                 */
                unsigned getHorizon() const;

                /**
                 * @brief This function will return the current set default value function.
                 *
                 * @return The currently set default value function.
                 */
                const ValueFunction & getValueFunction() const;

            private:
                // Parameters
                double discount_, epsilon_;
                unsigned horizon_;
                ValueFunction vParameter_;

                // Internals
                ValueFunction v1_;
                size_t S, A;

                // Internal methods
                /**
                 * @brief This function computes the single PRType of the MDP once for improved speed.
                 *
                 * @tparam M The type of the solvable MDP.
                 * @param m The MDP that needs to be solved.
                 *
                 * @return The Models's PRType.
                 */
                template <typename M, typename std::enable_if<is_model<M>::value, int>::type = 0>
                Table2D computeImmediateRewards(const M & model) const;

                /**
                 * @brief This function creates the Model's most up-to-date QFunction.
                 *
                 * @tparam M The type of the solvable MDP.
                 *
                 * @param m The MDP that needs to be solved.
                 * @param ir The immediate rewards of the model.
                 *
                 * @return A new QFunction.
                 */
                template <typename M, typename std::enable_if<is_model<M>::value, int>::type = 0>
                QFunction computeQFunction(const M & model, const Table2D & ir) const;

                /**
                 * @brief This function applies a single pass Bellman operator, improving the current ValueFunction estimate.
                 *
                 * This function computes the optimal value and action for
                 * each state, given the precomputed QFunction.
                 *
                 * @param q The precomputed QFunction.
                 * @param vOut The newly estimated ValueFunction.
                 */
                inline void bellmanOperator(const QFunction & q, ValueFunction * vOut) const;
        };

        template <typename M, typename std::enable_if<is_model<M>::value, int>::type>
        std::tuple<bool, ValueFunction, QFunction> ValueIteration::operator()(const M & model) {
            // Extract necessary knowledge from model so we don't have to pass it around
            S = model.getS();
            A = model.getA();
            discount_ = model.getDiscount();

            {
                // Verify that parameter value function is compatible.
                size_t size = std::get<VALUES>(vParameter_).size();
                if ( size != S ) {
                    if ( size != 0 )
                        std::cerr << "AIToolbox: Size of starting value function in ValueIteration::solve() is incorrect, ignoring...\n";
                    // Defaulting
                    v1_ = makeValueFunction(S);
                }
                else
                    v1_ = vParameter_;
            }

            auto ir = computeImmediateRewards(model);

            unsigned timestep = 0;
            double variation = epsilon_ * 2; // Make it bigger

            Values val0;
            QFunction q = makeQFunction(S, A);

            bool useEpsilon = checkDifferent(epsilon_, 0.0);
            while ( timestep < horizon_ && (!useEpsilon || variation > epsilon_) ) {
                ++timestep;

                auto & val1 = std::get<VALUES>(v1_);
                val0 = val1;

                q = computeQFunction(model, ir);
                bellmanOperator(q, &v1_);

                // We do this only if the epsilon specified is positive, otherwise we
                // continue for all the timesteps.
                if ( useEpsilon ) {
                    auto computeVariation = [](double lhs, double rhs) { return std::fabs(lhs - rhs); };
                    // We compute the difference and store it into v0 for comparison.
                    std::transform(std::begin(val1), std::end(val1), std::begin(val0), std::begin(val0), computeVariation);

                    variation = *std::max_element(std::begin(val0), std::end(val0));
                }
            }

            // We do not guarantee that the Value/QFunctions are the perfect ones, as we stop as within epsilon.
            return std::make_tuple(variation <= epsilon_, v1_, q);
        }

        template <typename M, typename std::enable_if<is_model<M>::value, int>::type>
        Table2D ValueIteration::computeImmediateRewards(const M & model) const {
            Table2D pr(boost::extents[S][A]);

            for ( size_t s = 0; s < S; ++s )
                for ( size_t a = 0; a < A; ++a )
                    for ( size_t s1 = 0; s1 < S; ++s1 )
                        pr[s][a] += model.getTransitionProbability(s,a,s1) * model.getExpectedReward(s,a,s1);

            return pr;
        }

        template <typename M, typename std::enable_if<is_model<M>::value, int>::type>
        QFunction ValueIteration::computeQFunction(const M & model, const Table2D & ir) const {
            QFunction q = ir;

            for ( size_t s = 0; s < S; ++s )
                for ( size_t a = 0; a < A; ++a )
                    for ( size_t s1 = 0; s1 < S; ++s1 )
                        q[s][a] += model.getTransitionProbability(s,a,s1) * discount_ * std::get<VALUES>(v1_)[s1];
            return q;
        }

        void ValueIteration::bellmanOperator(const QFunction & q, ValueFunction * v) const {
            auto & values  = std::get<VALUES> (*v);
            auto & actions = std::get<ACTIONS>(*v);

            for ( size_t s = 0; s < S; ++s ) {
                // Accessing an element like this creates a temporary. Thus we need to bind it.
                QFunction::const_reference ref = q[s];
                auto begin = std::begin(ref);
                auto it = std::max_element(begin, std::end(ref));

                values[s] = *it;
                actions[s] = std::distance(begin, it);
            }
        }
    }
}

#endif

Source file:

#include <AIToolbox/MDP/Algorithms/ValueIteration.hpp>

namespace AIToolbox {
    namespace MDP {
        ValueIteration::ValueIteration(unsigned horizon, double epsilon, ValueFunction v) : horizon_(horizon), vParameter_(v),
                                                                                            S(0), A(0)
        {
            setEpsilon(epsilon);
        }

        void ValueIteration::setEpsilon(double e) {
            if ( e < 0.0 ) throw std::invalid_argument("Epsilon must be >= 0");
            epsilon_ = e;
        }
        void ValueIteration::setHorizon(unsigned h) {
            horizon_ = h;
        }
        void ValueIteration::setValueFunction(ValueFunction v) {
            vParameter_ = v;
        }

        double                  ValueIteration::getEpsilon() const {
            return epsilon_;
        }
        unsigned                ValueIteration::getHorizon() const {
            return horizon_;
        }
        const ValueFunction &   ValueIteration::getValueFunction() const {
            return vParameter_;
        }
    }
}

Answer 1

I personally don't see the point in splitting things into hpp/cpp files here. For such a small amount of functionality that is basically just getters and setters, you may as well just place them in the header file, along with everything else.

---

I dislike the use of std::enable_if to toggle functions on and off when it isn't absolutely necessary. It does have a place, but this is generally when you need to perform overload resolution between multiple functions based on template parameters. If you simply want to ensure that some condition is met for a given template parameter, prefer to use static_assert. Without seeing exactly what the implementation of is_model is, it's hard to say exactly what this should be, but I imagine something like:

template <typename M>
QFunction ValueIteration::computeQFunction(const M & model, const Table2D & ir) const {
    static_assert(is_model<M>::value, "M must be a model!");
    //Implementation
    .....
}

On a separate note, you're missing an include here: there should be an #include <type_traits> for std::enable_if. You're also missing an #include <algorithm> (for things like std::max_element and std::transform), and an #include <cmath> for std::fabs. Even if other internal headers of yours are pulling these in, it's bad practice to rely on that - each file should pull in everything it needs by itself.

---

Perhaps the variable names S and A make sense in your particular context, but if possible, I'd try and give them more descriptive names.

---

Is there a reason that you've created an inner scope in ValueIterator::operator()?

template <typename M, typename std::enable_if<is_model<M>::value, int>::type>
std::tuple<bool, ValueFunction, QFunction> ValueIteration::operator()(const M & model) {
// Extract necessary knowledge from model so we don't have to pass it around
S = model.getS();
A = model.getA();
discount_ = model.getDiscount();

// Inner scope here?
{
    // Verify that parameter value function is compatible.
    size_t size = std::get<VALUES>(vParameter_).size();
    if ( size != S ) {
        if ( size != 0 )
            std::cerr << "AIToolbox: Size of starting value function in ValueIteration::solve() is incorrect, ignoring...\n";
        // Defaulting
        v1_ = makeValueFunction(S);
    }
    else
        v1_ = vParameter_;
}

Generally when I see this, I'm looking for something that uses RAII to clean thing up (files, locks, threads, etc). I can't see any particular reason for it here, which is a bit confusing when reading it.

---

Any particular reason why ValueFunction is passed by pointer here?

void ValueIteration::bellmanOperator(const QFunction & q, ValueFunction * v) const;

There are no check for nullptr, so I'd expect it to be passed by reference.

---

Mostly your code looks pretty reasonable. Using {} around for/if/etc has been mentioned, and I'll second that.

Answer 2

AI is not my area of knowledge, so I can't comment on the algorithm. The overall structure of the code seems good to me, so I'll just comment on a few style tidbits you could improve/change.

---

I recommend you always add bracings to controls statements (if, for, while, ...) even if it is a single line statement. Example:

for ( size_t s = 0; s < S; ++s )
{
    for ( size_t a = 0; a < A; ++a )
    {
        for ( size_t s1 = 0; s1 < S; ++s1 )
        {
            q[s][a] += model.getTransitionProbability(s,a,s1) * discount_ * std::get<VALUES>(v1_)[s1];
        }
    }
}

This way you will shield your code from the accidental insertion of an indented line that might look line it belongs to the statement:

for ( size_t s = 0; s < S; ++s )
    for ( size_t a = 0; a < A; ++a )
        for ( size_t s1 = 0; s1 < S; ++s1 )
            q[s][a] += model.getTransitionProbability(s,a,s1) * discount_ * std::get<VALUES>(v1_)[s1];
            if (somethingSomething(q))
            {
                ....
            }

On a quick pass through that block, you might not notice the mistake and think that the if will be executed inside the last for loop.

---

This line is quite long:

ValueIteration::ValueIteration(unsigned horizon, double epsilon, ValueFunction v) : horizon_(horizon), vParameter_(v),
                                                                                        S(0), A(0)
{
    setEpsilon(epsilon);
}

If you break the initialization list into other lines, it will be more readable:

ValueIteration::ValueIteration(unsigned int horizon, double epsilon, ValueFunction v) 
    : horizon_(horizon)
    , vParameter_(v)
    , S(0)
    , A(0)
{
    setEpsilon(epsilon);
}

---

The way you've indented these last lines (after the return type):

double                  ValueIteration::getEpsilon() const {
    return epsilon_;
}
unsigned                ValueIteration::getHorizon() const {
    return horizon_;
}
const ValueFunction &   ValueIteration::getValueFunction() const {
    return vParameter_;
}

That extra space doesn't seem to improve readability, I would just declare the functions normally with a single space after the return type:

double ValueIteration::getEpsilon() const {
    return epsilon_;
}

unsigned int ValueIteration::getHorizon() const {
    return horizon_;
}

const ValueFunction & ValueIteration::getValueFunction() const {
    return vParameter_;
}

---

Using unsigned by itself, like you do, seems a bit C-ish to me. I like to avoid the "implicit int" thing, so I recommend that you always use unsigned int instead.