I have this algorithm, called Variable Elimination, which is used quite often in my code-base (library). It is an algorithm for factor graphs (a type of bipartite graphs), used to efficiently maximize something over the nodes, thus I need to make it as fast as possible.

Since I had 5 separate implementations of this algorithm, and I'll probably need more in the future, I wanted to coalesce the algorithm into a single functor class.

Now, the main problem is that this is kind of a meta-algorithm. While the main steps are always the same, the specifics differ depending on what you are trying to actually accomplish by using the algorithm, and on what data structures you are using. Think of this as std::sort, which takes a lambda as a parameter for custom sorting; but here it's a lot messier. There's 4 nested loops total, and I may need to do custom operations between them. These custom operations may also need to store temporaries.

I have no idea whether there are best practices for this use-case, so the only idea I could come up with was to force the user to define a custom "global" structure which they have to pass in. This structure should provide a bunch of callbacks, which are then called within the algorithm itself to make everything run correctly. Some optional callbacks are called only if they exist using if constexpr.

I've already tested this, and while the usage is not super-pretty, I can use this algorithm in ~50 lines of code rather than 150 (and avoid duplication), which I guess it's a plus.

I'm looking for feedback on whether this approach makes any sense at all, and possibly what the alternatives are.


#include <AIToolbox/Factored/Utils/Core.hpp>
#include <AIToolbox/Factored/Utils/FactorGraph.hpp>

namespace AIToolbox::Factored {
     * @brief This class represents the Variable Elimination algorithm.
     * This class applies Variable Elimination to an input FactorGraph<Factor>.
     * Since the cross-sum steps in the algorithm differ from the type of node
     * in the graph, we require as input a separate structure which may contain
     * certain methods depending on what the use-case requires, and which holds
     * any needed temporaries to store for the duration of the algorithm.
     * In particular, this structure (the `global` parameter), *MUST* provide:
     * - A member `Factor newFactor` which stores the results of the cross-sum
     *   of each removed variable. At each iteration over the values of that
     *   variable's neighbors, we move from it, so be sure to re-initialize it
     *   if needed.
     * - A member `void crossSum(const Factor &)` function, which should
     *   perform the cross-sum of the input into the `newFactor` member
     *   variable.
     * - A member `void makeResult(FinalFactors &&)` method, which should
     *   process the final factors of the VE process in order to create your
     *   result.
     * Since VE usually requires custom computations, you can *OPTIONALLY*
     * define the following methods:
     * - A member `void beginRemoval()` method, which is called at the
     *   beginning of the removal of each variable.
     * - A member `void initNewFactor()` method, which is called when the
     *   `newFactor` variable needs to be initialized.
     * - A member `void beginCrossSum()` method, which is called at the
     *   beginning of each set of cross-sum operations.
     * - A member `void endCrossSum()` method, which is called at the end of
     *   each set of cross-sum operations.
     * - A member `bool isValidNewFactor()` method, which returns whether the
     *   `newFactor` variable can be used after all cross-sum operations.
     * - A member `void mergeRules(Rules &&, Rules &&)` method, which can be
     *   used to specify a custom step during the merge of the rules created by
     *   eliminating a variable with the previous ones.
     * All these functions can optionally be `const`; nothing changes. The
     * class will fail to compile if you provide a method with the required
     * name but with the wrong signature, as we would just skip it silently
     * otherwise.
     * @tparam Factor The Factor type to use.
    template <typename Factor>
    class GenericVariableElimination {
            using Rule = std::pair<PartialValues, Factor>;
            using Rules = std::vector<Rule>;
            using Graph = FactorGraph<Rules>;
            using FinalFactors = std::vector<Factor>;

             * @brief This operator performs the Variable Elimination operation on the inputs.
             * @param F The space of all factors to eliminate.
             * @param graph The already populated graph to perform VE onto.
             * @param global The global callback structure.
            template <typename Global>
            void operator()(const Factors & F, Graph & graph, Global & global);

             * @brief An helper struct to validate the interface of the global callback structure.
             * @tparam M The type of the global callback structure to validate.
            template <typename M>
            struct global_interface;

             * @brief This function removes the input factor from the graph.
             * @param F The space of all factors to eliminate.
             * @param graph The already populated graph to perform VE onto.
             * @param f The factor to eliminate.
             * @param finalFactors The storage of all the eliminated factors with no remaining neighbors.
             * @param global The global callback structure.
            template <typename Global>
            void removeFactor(const Factors & F, Graph & graph, const size_t f, FinalFactors & finalFactors, Global & global);

    template <typename Factor>
    template <typename M>
    struct GenericVariableElimination<Factor>::global_interface {
            #define STR2(X) #X
            #define STR(X) STR2(X)
            #define ARG(...) __VA_ARGS__

            // For each function we want to check, we are going to try each
            // overload in succession (char->int->long->...).
            // The first two simply accept the function with the approved
            // signature, whether it is const or not. The third checks whether
            // the member just exists, and reports that it probably has the
            // wrong signature (since we didn't match before).
            // The last just fails to find the match.
            #define MEMBER_CHECK(name, retval, input)                                       \
            template <typename Z> static constexpr auto name##Check(char) -> decltype(      \
                    static_cast<retval (Z::*)(input)> (&Z::name),                           \
                    bool()                                                                  \
                    ) { return true; }                                                      \
            template <typename Z> static constexpr auto name##Check(int) -> decltype(       \
                    static_cast<retval (Z::*)(input) const> (&Z::name),                     \
                    bool()                                                                  \
                    ) { return true; }                                                      \
            template <typename Z> static constexpr auto name##Check(long) -> decltype(      \
                    &Z::name,                                                               \
                    bool())                                                                 \
                    {                                                                       \
                        static_assert(!std::is_same_v<M, M>, "You provide a member '" STR(name) "' but with the wrong signature."); \
                        return false;                                                       \
                    }                                                                       \
            template <typename Z> static constexpr auto name##Check(...) -> bool { return false; }

            MEMBER_CHECK(beginRemoval, void, ARG(const Graph &, const typename Graph::FactorItList &, const typename Graph::VariableList &, const size_t))
            MEMBER_CHECK(initNewFactor, void, void)
            MEMBER_CHECK(beginCrossSum, void, void)
            MEMBER_CHECK(crossSum, void, const Factor &)
            MEMBER_CHECK(endCrossSum, void, void)
            MEMBER_CHECK(isValidNewFactor, bool, void)
            MEMBER_CHECK(mergeRules, Rules, ARG(Rules &&, Rules &&))
            MEMBER_CHECK(makeResult, void, FinalFactors &&)

            #undef MEMBER_CHECK
            #undef ARG
            #undef STR
            #undef STR2

            // All results are stored here for use later. All optional members
            // that do not exist, we simply do not call.
            enum {
                beginRemoval     = beginRemovalCheck<M>     ('\0'),
                initNewFactor    = initNewFactorCheck<M>    ('\0'),
                beginCrossSum    = beginCrossSumCheck<M>    ('\0'),
                crossSum         = crossSumCheck<M>         ('\0'),
                endCrossSum      = endCrossSumCheck<M>      ('\0'),
                isValidNewFactor = isValidNewFactorCheck<M> ('\0'),
                mergeRules       = mergeRulesCheck<M>       ('\0'),
                makeResult       = makeResultCheck<M>       ('\0'),

    template <typename Factor>
    template <typename Global>
    void GenericVariableElimination<Factor>::operator()(const Factors & F, Graph & graph, Global & global) {
        static_assert(global_interface<Global>::crossSum, "You must provide a crossSum method!");
        static_assert(global_interface<Global>::makeResult, "You must provide a makeResult method!");
        static_assert(std::is_same_v<Factor, decltype(global.newFactor)>, "You must provide a public 'Factor newFactor;' member!");

        FinalFactors finalFactors;

        // This can possibly be improved with some heuristic ordering
        while (graph.variableSize())
            removeFactor(F, graph, graph.variableSize() - 1, finalFactors, global);


    template <typename Factor>
    template <typename Global>
    void GenericVariableElimination<Factor>::removeFactor(const Factors & F, Graph & graph, const size_t f, FinalFactors & finalFactors, Global & global) {
        const auto factors = graph.getNeighbors(f);
        auto variables = graph.getNeighbors(factors);

        PartialFactorsEnumerator jointActions(F, variables, f);
        const auto id = jointActions.getFactorToSkipId();

        Rules newRules;

        if constexpr(global_interface<Global>::beginRemoval)
            global.beginRemoval(graph, factors, variables, f);

        // We'll now create new rules that represent the elimination of the
        // input variable for this round.
        const bool isFinalFactor = variables.size() == 1;

        while (jointActions.isValid()) {
            auto & jointAction = *jointActions;

            if constexpr(global_interface<Global>::initNewFactor)

            for (size_t sAction = 0; sAction < F[f]; ++sAction) {
                if constexpr(global_interface<Global>::beginCrossSum)

                jointAction.second[id] = sAction;
                for (const auto factor : factors)
                    for (const auto rule : factor->getData())
                        if (match(factor->getVariables(), rule.first, jointAction.first, jointAction.second))

                if constexpr(global_interface<Global>::endCrossSum)

            bool isValidNewFactor = true;
            if constexpr(global_interface<Global>::isValidNewFactor)
                isValidNewFactor = global.isValidNewFactor();

            if (isValidNewFactor) {
                if (!isFinalFactor) {
                    newRules.emplace_back(jointAction.second, std::move(global.newFactor));
                    // Remove new agent ID
                    newRules.back().first.erase(newRules.back().first.begin() + id);

        // And finally as usual in variable elimination remove the variable
        // from the graph and insert the newly created variable in.

        for (const auto & it : factors)

        if (!isFinalFactor && newRules.size()) {
            variables.erase(std::remove(std::begin(variables), std::end(variables), f), std::end(variables));

            auto newFactorNode = graph.getFactor(variables)->getData();

            if constexpr(global_interface<Global>::mergeRules) {
                newFactorNode = global.mergeRules(std::move(newFactorNode), std::move(newRules));
            } else {


Here is an example of a Global implementation, just to give you an idea:

struct Global {
    LP & lp;
    const size_t phiId;

    Factor newFactor;

    void initNewFactor() {
        newFactor = lp.row.size();

    void beginCrossSum() {
        lp.row[newFactor] = -1.0;
    void crossSum(const Factor & f) {
        lp.row[f] = 1.0;
    void endCrossSum() {
        lp.pushRow(LP::Constraint::LessEqual, 0.0);

        // Now do the reverse for all opposite rules (same rules +1)
        for (int i = lp.row.size() - 2; i >= 0; --i) {
            if (lp.row[i] != 0.0) {
                lp.row[i+1] = lp.row[i];
                lp.row[i] = 0.0;
        lp.pushRow(LP::Constraint::LessEqual, 0.0);
    void makeResult(VE::FinalFactors && finalFactors) {
        // Finally, add the two phi rules for all remaining factors.
        lp.row[phiId] = -1.0;

        for (const auto ruleId : finalFactors)
            lp.row[ruleId] = 1.0;

        lp.pushRow(LP::Constraint::LessEqual, 0.0);

        // Now do the reverse for all opposite rules (same rules +1)
        for (const auto ruleId : finalFactors) {
            lp.row[ruleId] = 0.0;
            lp.row[ruleId+1] = 1.0;

        lp.pushRow(LP::Constraint::LessEqual, 0.0);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.