graphnode.h

/* Spell-QTL  Software suite for the QTL analysis of modern datasets.
 * Copyright (C) 2016,2017  Damien Leroux <damien.leroux@inra.fr>, Sylvain Jasson <sylvain.jasson@inra.fr>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef _SPEL_BAYES_GRAPH_NODE_H_
#define _SPEL_BAYES_GRAPH_NODE_H_

#include "graphnode_base.h"

#include "../pedigree.h"

struct graph_type {
    node_vec rank;
    node_vec represented_by;
    std::vector<node_type> type;
    std::vector<colour_proxy> colour;
    std::vector<node_vec> neighbours_in;
    std::vector<node_vec> neighbours_out;
    std::vector<node_vec> inner_nodes;
    std::vector<var_vec> rules;
    var_vec variables;
    std::vector<var_vec> node_variables;
    std::map<variable_index_type, VariableIO> io;

    std::map<variable_index_type, node_index_type> interface_to_node;
    std::map<node_index_type, variable_index_type> node_to_interface;

    std::vector<std::shared_ptr<graph_type>> subgraphs;

    std::vector<std::shared_ptr<message_type>> tables;
    /*std::vector<message_type> state;*/
    std::map<var_vec, genotype_comb_type> domains;
    /* TODO suppress joint_parent_domains and all afferent code */
    std::map<var_vec, genotype_comb_type> joint_parent_domains;
    std::map<variable_index_type, char> ancestor_letters;
    /*std::vector<compute_state_operation_type> compute_state_ops;*/

    const graph_type* parent;
    node_index_type index_in_parent;

    bool aggregate_cycles;
    bool generate_interfaces;

    size_t n_alleles;

    /* THIS DOESN'T HAVE  TO BE SAVED/LOADED. TEMPORARY STATE IN ORDER TO COMPUTE THE SEQUENCES OF OPERATIONS. */
    std::vector<var_vec> annotations;

    std::vector<bool> is_dh;

    graph_type& operator = (graph_type&& other) = delete;

    /* node index 0 is the trash bin */

    graph_type(const graph_type& other) = delete;

    graph_type()
        : rank(1), represented_by(1), type(1), colour(1), neighbours_in(1), neighbours_out(1), inner_nodes(1), rules(1), variables(1), node_variables(1), io(), interface_to_node(), node_to_interface(), subgraphs(1), tables(1), /*state(1),*/ parent(nullptr), index_in_parent(0),aggregate_cycles(true), generate_interfaces(true), n_alleles(1), annotations(1), is_dh(1, false)
    {}

    graph_type(size_t n_al)
        : rank(1), represented_by(1), type(1), colour(1), neighbours_in(1), neighbours_out(1), inner_nodes(1), rules(1), variables(1), node_variables(1), io(), interface_to_node(), node_to_interface(), subgraphs(1), tables(1), /*state(1),*/ parent(nullptr), index_in_parent(0), aggregate_cycles(true), generate_interfaces(true), n_alleles(n_al), annotations(1), is_dh(1, false)
    {}

    bool is_aggregate(node_index_type node) const { return type[node] == Aggregate/*inner_nodes[node].size() > 1*/; }

    bool
        is_compound_interface(node_index_type node) const
        {
            if (is_aggregate(node)) {
                for (node_index_type n: inner_nodes[node]) {
                    if (type[n] != Interface) {
                        return false;
                    }
                }
                return true;
            }
            return false;
        }

    bool
        is_interface(node_index_type node) const
        {
            if (is_aggregate(node)) {
                for (node_index_type n: inner_nodes[node]) {
                    if (type[n] != Interface) {
                        return false;
                    }
                }
                return true;
            }
            return type[node] == Interface;
        }

    bool
        is_computable(node_index_type node) const
        {
            return !is_interface(node);
        }

    size_t size() const { return rank.size(); }

    node_vec
        active_nodes() const
        {
            node_vec ret;
            ret.reserve(represented_by.size());
            for (node_index_type i = 1; i < represented_by.size(); ++i) {
                if (represented_by[i] == i) {
                    ret.push_back(i);
                }
            }
            return ret;
        }

    node_vec
        resolve_vector(const node_vec& vec) const
        {
            node_vec ret;
            ret.reserve(vec.size());
            for (node_index_type i: vec) {
                node_index_type r = resolve(i);
                if (r > 0) {
                    ret.push_back(r);
                }
            }
            sort_and_unique(ret);
            return ret;
        }

    var_vec
        interface_nodes(var_vec inputs) const
        {
            var_vec ret;
            for (variable_index_type v: inputs) {
                ret.push_back(resolve(interface_to_node.find(v)->second));
            }
            sort_and_unique(ret);
            return ret;
        }

    std::vector<edge_type>
        active_edges() const
        {
            std::vector<edge_type> ret;
            for (node_index_type n: active_nodes()) {
                for (node_index_type o: nei_out(n)) {
                    ret.emplace_back(this, n, o);
                }
            }
            return ret;
        }

    node_vec
        nei_in(node_index_type n) const
        {
            return resolve_vector(neighbours_in[n]);
        }

    node_vec
        nei_out(node_index_type n) const
        {
            return resolve_vector(neighbours_out[n]);
        }

    node_vec
        all_nei(node_index_type n) const
        {
            return nei_in(n) + nei_out(n);
        }

    void
        remove_link(node_index_type nin, node_index_type nout)
        {
            neighbours_out[nin] = resolve_vector(neighbours_out[nin]) - node_vec{nout};
            neighbours_in[nout] = resolve_vector(neighbours_in[nout]) - node_vec{nin};
        }

    void
        dump_node(node_index_type n)
        {
            MSG_DEBUG(
                   '[' << rank[n] << "] " << (is_interface(n) ? "INTERFACE " : (n == inner_nodes[n][0] ? "FACTOR " : "AGGREGATE ")) << n << std::endl
                << "  creation rule " << rules[n] << std::endl
                << "  represented by " << represented_by[n] << " (" << resolve(n) << ')' << std::endl
                << "  colour " << get_colour_impl(colour[n]) << std::endl
                << "  inputs " << neighbours_in[n] << " (" << resolve_vector(neighbours_in[n]) << ')' << std::endl
                << "  outputs " << neighbours_out[n] << " (" << resolve_vector(neighbours_out[n]) << ')' << std::endl
                << "  inner nodes " << inner_nodes[n] << std::endl
                << "  variable(s) " << variables_of(n) << std::endl
                );
            if (node_domains.size() == rank.size()) {
                MSG_DEBUG("  TABLE " << node_domains[n] << std::endl);
            }
            if (inner_nodes[n].size() > 1) {
                if (subgraphs[n]) {
                    scoped_indent _("  | ");
                    subgraphs[n]->dump_active();
                }
            }
            MSG_DEBUG("");
        }

#define DUMP_SZ(_x) << " * " #_x " " << _x.size() << std::endl

    void
        dump_sizes() const
        {
            MSG_DEBUG(""
                DUMP_SZ(rank)
                DUMP_SZ(type)
                DUMP_SZ(rules)
                DUMP_SZ(colour)
                DUMP_SZ(variables)
                DUMP_SZ(inner_nodes)
                DUMP_SZ(neighbours_in)
                DUMP_SZ(neighbours_out)
                DUMP_SZ(represented_by)
                );
        }

    void
        dump()
        {
            MSG_DEBUG("ALL NODES");
            /*dump_sizes();*/
            for (node_index_type i = 1; i < rank.size(); ++i) {
                dump_node(i);
            }
        }

    void
        dump_active()
        {
            MSG_DEBUG("ACTIVE NODES");
            /*dump_sizes();*/
            for (node_index_type i: active_nodes()) {
                dump_node(i);
            }
        }

    void
        compute_ranks()
        {
            std::vector<bool> visited(rank.size(), false);
            compute_ranks(active_nodes(), visited);
        }

    void
        compute_ranks(const node_vec& nodes, std::vector<bool>& visited)
        {
            for (node_index_type n: nodes) {
                if (visited[n]) { continue; }
                visited[n] = true;
                auto nin = nei_in(n);
                compute_ranks(nin, visited);
                rank[n] = 0;
                if (nin.size()) {
                    for (node_index_type i: nin) {
                        if (rank[n] < rank[i]) {
                            rank[n] = rank[n];
                        }
                    }
                    ++rank[n];
                }
            }
        }

    typedef std::pair<size_t, var_vec> emitter_and_interface_type;
    struct compare_eai {
        bool operator () (const emitter_and_interface_type& e1, const emitter_and_interface_type& e2) const { return e1.first < e2.first || (e1.first == e2.first && e1.second < e2.second); }
    };

    typedef std::map<emitter_and_interface_type, size_t, compare_eai> interface_map_type;

    node_index_type
        create_interface(const var_vec& varset)
        {
            node_index_type i;
            if (varset.size() == 1) {
                i = add_interface(node_vec{}, varset.front());
                colour[i] = create_colour();
            } else {
                node_vec iv;
                iv.reserve(varset.size());
                for (variable_index_type v: varset) {
                    iv.push_back(add_interface(node_vec{}, v));
                }
                i = add_node(node_vec{}, node_vec{}, var_vec{}, create_colour(), Aggregate, iv, -1);
                for (node_index_type ni: iv) {
                    represented_by[ni] = i;
                    colour[ni] = colour[i];
                }
            }
            return i;
        }

    node_index_type
        create_interface_between(node_index_type n1, node_index_type n2)
        {
            auto varset = variables_of(n1) % variables_of(n2);
            node_index_type i = create_interface(varset);
            colour[i] = get_colour_impl(colour[n1]);
            neighbours_in[i].push_back(n1);
            neighbours_out[i].push_back(n2);
            return i;
        }

    void
        rebuild_interface_between(node_index_type n1, node_index_type n2, interface_map_type& interface_map)
        {
            auto varset = variables_of(n1) % variables_of(n2);
            node_index_type& i = interface_map[{n1, varset}];
            if (i == 0) {
                i = create_interface_between(n1, n2);
            } else {
                neighbours_out[i] = neighbours_out[i] + node_vec{n2};
            }
            neighbours_in[n2] = nei_in(n2) + node_vec{i} - nei_in(i);
            neighbours_out[n1] = nei_out(n1) + node_vec{i} - nei_out(i);
            remove_link(n1, n2);
            /*filter_out_and_replace_by(neighbours_out[n1], node_vec{n2}, i);*/
            /*filter_out_and_replace_by(neighbours_in[n2], node_vec{n1}, i);*/
            /*MSG_DEBUG("REBUILD " << n1 << ' ' << i << ' ' << n2);*/
            /*dump_node(n1);*/
            /*dump_node(i);*/
            /*dump_node(n2);*/
        }

    void
        remove_nei(node_vec& neighbours, node_index_type n)
        {
            node_vec new_nei;
            new_nei.reserve(neighbours.size());
            auto N = resolve(n);
            for (auto x: resolve_vector(neighbours)) {
                if (x != N) {
                    new_nei.push_back(x);
                }
            }
            new_nei.swap(neighbours);
        }

    bool
        suppress_all_cycles()
        {
            bool redo;
            bool redone = false;
            do {
                redo = false;
                for (node_index_type n: active_nodes()) {
                    auto nei = nei_in(n);
                    /*MSG_DEBUG("LOOK FOR CYCLES FROM NEIGHBOURS OF " << n << "   " << nei);*/
                    for (auto i = nei.begin(), j = nei.end(); i != j; ++i) {
                        std::list<node_index_type> path;
                        for (auto k = i + 1; k != j; ++k) {
                            path = find_path_between_parents(*i, *k, n);
                            if (path.size()) {
                                /*MSG_DEBUG("FOUND BETWEEN " << (*i) << " AND " << (*k) << ": " << path);*/
                                /*for (node_index_type x: path) {*/
                                    /*dump_node(x);*/
                                /*}*/
                                break;
                            }
                        }
                        if (path.size()) {
                            MSG_DEBUG("FOUND CYCLE ENDING ON " << n << " === " << path);
                            path.push_back(n);
                            MSG_DEBUG("path to aggregate " << path);
                            aggregate_path(path);
                            redo = true;
                            redone = true;
                            break;
                        }
                    }
                }
            } while (redo);
            return redone;
        }

    void
        finalize_stage1(bool reconstruct_itf=true)
        {
            /*check_neighbours();*/
            debug_graph("finalize", 10);
            /*scoped_indent _("[stage1] ");*/
            /* search for cycles and aggregate them */
            interface_map_type interface_map;
            if (reconstruct_itf) {
                /*dump_active();*/
                auto edges = active_edges();
                for (const auto& e: edges) {
                    if (is_interface(e.first) ^ !is_interface(e.second)) {
                        /* By construction, it's a factor->factor edge, never an interface->interface edge */
                        rebuild_interface_between(e.first, e.second, interface_map);
                    }
                }
            }
            /*check_neighbours();*/
            debug_graph("finalize", 11);
            suppress_all_cycles();
            /*check_neighbours();*/
            debug_graph("finalize", 12);

            if (size()) {
                node_vec A = active_nodes();

                for (node_index_type n: A) {
                    if (is_aggregate(n) && !is_interface(n)) {
                        /*MSG_DEBUG("subgraphs.size() = " << subgraphs.size() << " n = " << n);*/
                        subgraphs[n] = subgraph(n);
                        /* protect all interfaced variables with this level */
                        var_vec var;
                        for (node_index_type ni: nei_in(n)) {
                            auto nv = variables_of(ni);
                            var.insert(var.end(), nv.begin(), nv.end());
                        }
                        for (node_index_type no: nei_out(n)) {
                            auto nv = variables_of(no);
                            var.insert(var.end(), nv.begin(), nv.end());
                        }
                        sort_and_unique(var);
                        for (auto v: var) {
                            subgraphs[n]->io[v] = Input | Output;
                        }
                    }
                    /* then optimize */
                    /*subgraphs[n]->optimize();*/
                }
                debug_graph("finalize", 13);
            }
        }

    void
        finalize_stage2(bool reconstruct_itf)
        {
            if (reconstruct_itf) {
                /* Create interfaces for all variables used in this layer not already represented in a top-level interface */
                var_vec all_variables;
                std::vector<bool> var_represented(1 + *std::max_element(variables.begin(), variables.end()), false);
                node_vec A = active_nodes();
                for (node_index_type n: A) {
                    if (is_interface(n)) {
                        for (variable_index_type v: variables_of(n)) {
                            var_represented[v] = true;
                        }
                    } else if (is_aggregate(n) && nei_in(n).size() == 0) {
                        var_vec vv;
                        for (node_index_type subn: subgraphs[n]->active_nodes()) {
                            if (subgraphs[n]->is_interface(subn) && subgraphs[n]->nei_in(subn).size() == 0) {
                                vv = vv + subgraphs[n]->variables_of(subn);
                            }
                        }
                        if (vv.size()) {
                            node_index_type i = create_interface(vv);
                            neighbours_out[i].push_back(n);
                            neighbours_in[n] = nei_in(n) + node_vec{i};
                            colour[i] = get_colour_impl(colour[n]);
                            for (variable_index_type v: vv) {
                                var_represented[v] = true;
                            }
                        }
                    }
                }
                debug_graph("finalize", 21);

                /* Discover input interfaces in order to merge distinct trees */
                /* Actually we need to discover their neighbours so we can sort-unique them by colours */
                if (parent != NULL) {
                    node_vec itf_nei;
                    for (node_index_type n: A) {
                        if (is_interface(n) && nei_in(n).size() == 0) {
                            itf_nei = itf_nei + nei_out(n);
                        } else if (!is_interface(n)) {
                            for (variable_index_type v: variables_of(n)) {
                                if (var_represented[v]) { continue; }
                                node_index_type i = add_interface(node_vec{n}, v);
                                nei_out(n).push_back(i);
                            }
                        }
                    }
                    std::sort(itf_nei.begin(), itf_nei.end(), [this] (node_index_type a, node_index_type b) { return colour[a] < colour[b]; });
                    itf_nei.erase(std::unique(itf_nei.begin(), itf_nei.end(), [this] (node_index_type a, node_index_type b) { return colour_equal(colour[a], colour[b]); }), itf_nei.end());

                    node_vec itf;
                    node_vec out, inner;
                    for (node_index_type nei: itf_nei) {
                        for (node_index_type n: nei_in(nei)) {
                            if (nei_in(n).size() == 0) {
                                itf.push_back(n);
                                out = out + nei_out(n);
                                inner = inner + inner_nodes[n];
                            }
                        }
                    }

                    /* now also add single interfaces that have no neighbours */

                    for (node_index_type i: A) {
                        if (all_nei(i).size() == 0) {
                            itf.push_back(i);
                            inner = inner + inner_nodes[i];
                        }
                    }
                    std::sort(itf.begin(), itf.end());
                    std::sort(inner.begin(), inner.end());

                    /* Merge inputs of distinct trees to make one single tree */
#if 1
                    /*MSG_DEBUG("Input interfaces " << itf);*/
                    /*for (auto i: itf) {*/
                        /*MSG_DEBUG(" * " << i << " has colour " << get_colour_impl(colour[i]));*/
                    /*}*/
                    /*std::sort(itf.begin(), itf.end(), [this] (node_index_type a, node_index_type b) { return colour[a] < colour[b]; });*/
                    /*itf.erase(std::unique(itf.begin(), itf.end(), [this] (node_index_type a, node_index_type b) { return colour_equal(colour[a], colour[b]); }), itf.end());*/
                    /*MSG_DEBUG("Colour unique input interfaces " << itf);*/
                    /*for (auto i: itf) {*/
                        /*MSG_DEBUG(" * " << i << " has colour " << get_colour_impl(colour[i]));*/
                    /*}*/
                    if (itf.size() > 1) {
                        /* FIXME merge interface IIF factor colours are different. Otherwise this creates a cycle. */
                        node_index_type agr = add_node(node_vec{}, out, var_vec{}, create_colour(), Aggregate, inner, -1);
                        node_vec va = {agr};
                        for (node_index_type i: itf) {
                            represented_by[i] = agr;
                        }
                        for (node_index_type n: out) {
                            neighbours_in[n] = nei_in(n) - itf + va;
                            if (!colour_equal(colour[n], colour[agr])) {
                                assign_colour_impl(colour[n], colour[agr]);
                            }
                        }
                    }
                    debug_graph("finalize", 22);
#endif
                /*} else {*/
                    /* FIXME we shouldn't have to recompute this here */
                    A = active_nodes();
                    for (node_index_type n: A) {
                        if (is_interface(n)) {
                            for (variable_index_type v: variables_of(n)) {
                                var_represented[v] = true;
                            }
                        }
                    }
                    for (node_index_type n: A) {
                        if (!is_interface(n)) {
                            var_vec all_output;
                            for (node_index_type o: nei_out(n)) {
                                all_output = all_output + variables_of(o);
                            }
                            for (variable_index_type v: variables_of(n)) {
                                if (var_represented[v] || std::find(all_output.begin(), all_output.end(), v) == all_output.end()) { continue; }
                                node_index_type i = add_interface(node_vec{n}, v);
                                neighbours_out[n].push_back(i);
                                colour[i] = colour[n];
                            }
                        }
                    }
                    debug_graph("finalize", 23);
                }
            }
            /* unbox single-factor layers */
            for (node_index_type n: active_nodes()) {
                for (node_index_type o: nei_out(n)) {
                    if (!colour_equal(colour[n], colour[o])) {
                        assign_colour_impl(colour[o], colour[n]);
                    }
                }
            }
#if 1
            for (node_index_type n: active_nodes()) {
                if (type[n] == Aggregate && !is_interface(n)) {
                    size_t count_factors = 0;
                    node_index_type factor_node = (node_index_type) -1;
                    for (node_index_type i: inner_nodes[n]) {
                        if (type[i] == Factor) {
                            ++count_factors;
                            factor_node = i;
                        }
                    }
                    MSG_DEBUG("Aggregate #" << n << " has " << count_factors << " factors");
                    if (count_factors == 1) {
                        represented_by[factor_node] = factor_node;
                        neighbours_in[factor_node] = nei_in(n);
                        neighbours_out[factor_node] = nei_out(n);
                        represented_by[n] = factor_node;
                        subgraphs[n].reset();
                        for (node_index_type i: nei_in(n)) {
                            neighbours_out[i] = nei_out(i);
                        }
                        for (node_index_type o: nei_out(n)) {
                            neighbours_in[o] = nei_in(o);
                        }
                    }
                }
            }
#endif
        }

    void
        finalize_import_inputs()
        {
            for (node_index_type n: active_nodes()) {
                if (is_aggregate(n) && !is_interface(n)) {
                    std::vector<var_vec> inputs;
                    for (node_index_type i: nei_in(n)) {
                        inputs.push_back(variables_of(i));
                    }
                }
            }
        }

    void
        finalize(bool reconstruct_itf=true)
        {
            scoped_indent _("[finalize] ");
            finalize_stage1(reconstruct_itf);
            finalize_stage2(reconstruct_itf);
            update_all_ranks();
            /*if (reconstruct_itf) {*/
                /*finalize_import_inputs();*/
            /*}*/
        }

    node_index_type
        add_node(const node_vec& in, const node_vec& out,
                 const var_vec& rule, colour_proxy col, node_type t,
                 const node_vec& inner, variable_index_type var)
        {
            node_index_type ret = rank.size();
            /*MSG_DEBUG("adding node " << ret);*/
            neighbours_out.emplace_back(out);
            neighbours_in.emplace_back(in);
            rules.emplace_back(rule);
            colour.emplace_back(col);
            if (in.size()) {
                size_t r = 0;
                for (node_index_type i: in) {
                    r = std::max(r, rank[i]);
                }
                rank.push_back(r + 1);
            } else {
                rank.push_back(0);
            }
            /*MSG_DEBUG("rank=" << rank.back());*/
            type.push_back(t);
            variables.push_back(var);
            node_variables.emplace_back();
            represented_by.push_back(ret);
            subgraphs.emplace_back();
            tables.emplace_back();
            /*state.emplace_back();*/
            annotations.emplace_back();
            is_dh.push_back(false);
            if (inner.size()) {
                inner_nodes.emplace_back(inner);
            } else {
                inner_nodes.emplace_back(node_vec{ret});
            }
            /*dump();*/
            return ret;
        }

    node_index_type
        resolve_interface(variable_index_type var)
        {
            auto it = interface_to_node.find(var);
            if (it == interface_to_node.end()) {
                /*MSG_DEBUG("resolve_interface(" << var << ") => new interface");*/
                if (generate_interfaces) {
                    return add_interface(node_vec{}, var);
                } else {
                    return (node_index_type) -1;
                }
            }
            /*MSG_DEBUG("resolve_interface(" << var << ") => resolve(" << it->second << ") = " << resolve(it->second));*/
            return resolve(it->second);
        }

    node_index_type
        add_interface(const node_vec& producer, variable_index_type var, bool force=false)
        {
            node_index_type ret = add_node(producer, node_vec{}, var_vec{}, create_colour(), Interface, node_vec{}, var);
            if (force || interface_to_node.find(var) == interface_to_node.end()) {
                interface_to_node[var] = ret;
                node_to_interface[ret] = var;
            }
            for (node_index_type p: producer) {
                neighbours_out[p].push_back(ret);
            }
            return ret;
        }

    node_index_type
        add_factor(const var_vec& rule, colour_proxy col,
                 variable_index_type var)
        {
            node_vec in, out;
            /*
            if (generate_interfaces) {
                for (variable_index_type v: rule) {
                    in.push_back(resolve_interface(v));
                }
            } else {
                for (variable_index_type v: rule) {
                    auto it = interface_to_node.find(v);
                    if (it != interface_to_node.end()) {
                        in.push_back(resolve(it->second));
                    }
                }
            }
            sort_and_unique(in);
            */
            node_index_type ret = add_node(in, node_vec{}, rule, col, Factor, node_vec{}, var);
            for (node_index_type n: in) {
                neighbours_out[n].push_back(ret);
            }
            if (generate_interfaces) {
                node_index_type i = add_node(node_vec{ret}, node_vec{}, var_vec{}, colour[ret], Interface, node_vec{}, var);
                interface_to_node[var] = i;
                node_to_interface[i] = var;
                neighbours_out[ret].push_back(i);
            } else {
                interface_to_node[var] = ret;
                node_to_interface[ret] = var;
            }
            return ret;
        }

    node_index_type
        add_factor(variable_index_type var)
        {
            /*MSG_DEBUG("add_factor(" << var << ')');*/
            node_index_type ret = add_factor(var_vec{}, create_colour(), var);
            /*dump();*/
            return ret;
        }

    node_index_type
        resolve(node_index_type n) const
        {
            while (n != represented_by[n]) { n = represented_by[n]; }
            return n;
        }

    node_index_type
        add_factor(variable_index_type v1, variable_index_type var)
        {
            /*MSG_DEBUG("add_factor(" << v1 << ", " << var << ')');*/
            node_index_type p1r;
            if (!generate_interfaces) {
                auto it = interface_to_node.find(v1);
                if (it == interface_to_node.end()) {
                    auto ret = add_factor(var);
                    rules.back() = {v1};
                    return ret;
                } else {
                    p1r = resolve(it->second);
                }
            } else {
                p1r = resolve_interface(v1);
            }
            node_index_type ret = add_factor(var_vec{v1}, colour[p1r], var);
            neighbours_in[ret] = {p1r};
            neighbours_out[p1r] = {ret};
            compute_ranks();
            /*dump();*/
            return ret;
        }

    void
        aggregate_path(std::list<node_index_type>& path)
        {
            parents_of_max_rank aggr_first;
            while (path.size() > 2 && ((aggr_first = find_parents_of_max_rank(path)), !aggregate(*aggr_first.p1, *aggr_first.p2, *aggr_first.child, path, aggr_first.p1))) {
                path.erase(aggr_first.p1);
                path.erase(aggr_first.p2);
                path.erase(aggr_first.child);
            }
        }

    node_index_type
        add_factor(variable_index_type v1,variable_index_type v2, variable_index_type var)
        {
            /*MSG_DEBUG("add_factor(" << v1 << ", " << v2 << ", " << var << ')');*/
            node_index_type p1r, p2r;
            /* first, search for a factor/aggregate containing both parents */
            var_vec sorted_rule = v1 < v2 ? var_vec{v1, v2} : var_vec{v2, v1};
            bool found = false;
            node_index_type common = 0;
            for (auto n: active_nodes()) {
                if (type[n] != Factor) {
                    continue;
                }
                auto varset = variables_of(n);
                auto result = sorted_rule % varset;
                MSG_DEBUG("searching for rule " << sorted_rule << " in #" << n << ' ' << varset << " => " << result);
                if (result == sorted_rule) {
                    MSG_DEBUG("Found!");
                    found = true;
                    common = n;
                    break;
                }
            }
            if (!generate_interfaces) {
                if (found) {
                    p1r = p2r = common;
                } else {
                    auto it = interface_to_node.find(v1);
                    if (it == interface_to_node.end()) {
                        node_index_type ret = add_factor(v2, var);
                        rules.back() = {v1, v2};
                        return ret;
                    } else {
                        p1r = resolve(it->second);
                    }
                    it = interface_to_node.find(v2);
                    if (it == interface_to_node.end()) {
                        node_index_type ret = add_factor(v1, var);
                        rules.back() = {v1, v2};
                        return ret;
                    } else {
                        p2r = resolve(it->second);
                    }
                }
            } else {
                if (found) {
                    found = false;
                    node_index_type agr;
                    for (auto n: nei_out(common)) {
                        if (sorted_rule % variables_of(n) == sorted_rule) {
                            agr = n;
                            found = true;
                            break;
                        }
                    }
                    if (!found) {
                        node_vec common_fac = {common};
                        node_index_type i1, i2;
                        i1 = add_interface(common_fac, v1);
                        i2 = add_interface(common_fac, v2);
                        agr = add_node(common_fac, node_vec{}, var_vec{}, colour[common_fac.front()], Aggregate, node_vec{i1, i2}, -1);
                        represented_by[i1] = represented_by[i2] = agr;
                    }
                    p1r = p2r = agr;
                } else {
                    p1r = resolve_interface(v1);
                    p2r = resolve_interface(v2);
                }

                /*p1r = resolve_interface(v1);*/
                /*p2r = resolve_interface(v2);*/
            }

            MSG_DEBUG("p1r=" << p1r << " p2r=" << p2r);

            if (p1r > p2r) {
                p1r ^= p2r;
                p2r ^= p1r;
                p1r ^= p2r;
            }

            node_index_type ret;

            if (p1r != (node_index_type) -1) {
                ret = add_factor(var_vec{v1, v2}, colour[p1r], var);

                bool cycle = false;
                if (p1r != p2r) {
                    neighbours_in[ret] = {p1r, p2r};
                    neighbours_out[p1r].push_back(ret);
                    neighbours_out[p2r].push_back(ret);
                    if (colour_equal(colour[p1r], colour[p2r])) {
                        /* cycle! */
                        MSG_DEBUG("cycle!");
                        cycle = true;
                    } else {
                        MSG_DEBUG("no cycle.");
                        assign_colour_impl(colour[p1r], colour[p2r]);
                    }
                } else {
                    neighbours_in[ret] = {p1r};
                    neighbours_out[p1r].push_back(ret);
                }
                rank[ret] = std::max(rank[p1r], rank[p2r]) + 1;
                dump();
                MSG_QUEUE_FLUSH();

                if (cycle && aggregate_cycles) {
                    MSG_DEBUG("search a path between " << p1r << " and " << p2r);
                    std::vector<bool> visited(rank.size(), false);
                    visited[ret] = true;
                    auto path = find_shortest_path(node_vec{p1r, p2r},
                            [&, this](node_index_type n) { return nei_in(n) + nei_out(n); },
                            [&] (node_index_type n) { return visited[n]; },
                            [&] (node_index_type n) { visited[n] = true; });
                    path.push_back(ret);
                    MSG_DEBUG("New algo path: " << path);
                    path = find_path_between_parents(p1r, p2r, ret);
                    MSG_DEBUG("Old algo path: " << path);
                    /*dump_active();*/
                    MSG_DEBUG("Found path: "; for (size_t n: path) { std::cout << ' ' << n; } std::cout);
                    /*size_t min_rank = 1 + rank[*std::min_element(path.begin(), path.end(), [this](node_index_type i1, node_index_type i2) { return rank[i1] < rank[i2]; })];*/
                    aggregate_path(path);
                }
            } else {
                ret = add_factor(var_vec{v1, v2}, create_colour(), var);
            }

            compute_ranks();
            /*dump();*/
            dump_active();
            return ret;
        }

    void
        add_ancestor(variable_index_type id)
        {
            auto& domain = domains[{id}];
            char letter = ancestor_letters.size() + 'a';
            ancestor_letters[id] = letter;
            for (char al = 0; al < (char) n_alleles; ++al) {
                domain.m_combination.emplace_back(genotype_comb_type::element_type{{{id, {letter, letter, al, al}}}, 1.});
            }
            /*MSG_DEBUG("Domains " << domains);*/
        }

    node_index_type
        add_cross(variable_index_type p1, variable_index_type p2, variable_index_type id)
        {
            /*MSG_DEBUG("################ ADD FACTOR " << p1 << " " << p2 << " " << id);*/
            /*auto factor = compute_factor_table(id, var_vec{p1, p2}, false);*/
            node_index_type ret = add_factor(p1, p2, id);
            /*tables[ret] = factor;*/
            /*state[ret] = *tables[ret];*/
            /*MSG_DEBUG("Computed factor for #" << id << ": " << (*tables[ret]));*/
            /*operations.push_back(op);*/
            return ret;
        }

    node_index_type
        add_dh(variable_index_type p1, variable_index_type id)
        {
            /*auto factor = compute_factor_table(id, var_vec{p1}, true);*/
            node_index_type ret = add_factor(p1, id);
            /*tables[ret] = factor;*/
            /*state[ret] = *tables[ret];*/
            is_dh[ret] = true;
            /*MSG_DEBUG("Computed factor for #" << ret << ": " << (*tables[ret]));*/
            /*operations.push_back(op);*/
            return ret;
        }

    node_index_type
        add_selfing(variable_index_type p1, variable_index_type id)
        {
            /*auto factor = compute_factor_table(id, var_vec{p1}, false);*/
            node_index_type ret = add_factor(p1, id);
            /*tables[ret] = factor;*/
            /*state[ret] = *tables[ret];*/
            /*MSG_DEBUG("Computed factor for #" << ret << ": " << (*tables[ret]));*/
            /*operations.push_back(op);*/
            return ret;
        }

    bool
        find_ascending_path(node_index_type p1, node_index_type p2, node_vec& path, std::vector<bool>& visited)
        {
            /*MSG_DEBUG("path from " << p1 << " to " << p2);*/
            if (p1 == p2) {
                return true;
            }
            for (node_index_type i: nei_in(p1)) {
                i = resolve(i);
                if (visited[i]) {
                    return false;
                }
                visited[i] = true;
                if (find_ascending_path(i, p2, path, visited)) {
                    path.push_back(i);
                    return true;
                }
            }
            return false;
        }

    std::vector<bool>
        create_visited()
        {
            return std::vector<bool>(rank.size(), false);
        }

    node_vec
        find_aggregate_chain(node_index_type p1, node_index_type p2)
        {
            node_vec path;
            auto vis1 = create_visited(), vis2 = create_visited();

            if (rank[p1] > rank[p2] && find_ascending_path(p1, p2, path, vis1)) {
                /*MSG_DEBUG("found ascending path from " << p1 << " to " << p2);*/
                path.push_back(p1);
            } else if (rank[p2] > rank[p1] && find_ascending_path(p2, p1, path, vis2)) {
                /*MSG_DEBUG("found ascending path from " << p2 << " to " << p1);*/
                path.push_back(p2);
            } else {
                /*MSG_DEBUG("no ascending path between parents");*/
                path = {p1, p2};
            }
            return path;
        }

    bool
        filter_out(node_vec& vec, const node_vec& aggr)
        {
            node_vec tmp(vec.size());
            auto it = std::set_difference(vec.begin(), vec.end(), aggr.begin(), aggr.end(), tmp.begin());
            tmp.resize(it - tmp.begin());
            tmp.swap(vec);
            return tmp.size() != vec.size();
        }

    void
        filter_out_and_replace_by(node_vec& vec, const node_vec& aggr, node_index_type new_node)
        {
            /*MSG_DEBUG("filter_out_and_replace_by(" << vec << ", " << aggr << ", " << new_node << ')');*/
            /*if (filter_out(vec, aggr)) {*/
                /*vec.push_back(new_node);*/
            /*}*/
            node_vec tmp;
            tmp.reserve(vec.size());
            for (auto n: vec) {
                node_index_type nr = resolve(n);
                /*MSG_DEBUG("" << nr << " vs " << aggr);*/
                if (std::find(aggr.begin(), aggr.end(), nr) == aggr.end()) {
                    tmp.push_back(n);
                }
            }
            sort_and_unique(tmp);
            if (tmp.size() != vec.size()) {
                tmp.push_back(new_node);
            }
            vec.swap(tmp);
            /*MSG_DEBUG(" => " << vec);*/
        }

    bool
        update_rank(node_index_type node)
        {
            node_index_type r = 0;
            if (nei_in(node).size()) {
                for (node_index_type i: nei_in(node)) {
                    if (rank[i] > r) {
                        r = rank[i];
                    }
                }
                ++r;
            }
            if (rank[node] != r) {
                rank[node] = r;
                return true;
            }
            return false;
        }

    void
        propagate_update_rank(node_index_type source)
        {
            auto nout = nei_out(source);
            std::deque<node_index_type> stack(nout.begin(), nout.end());
            while (stack.size()) {
                node_index_type n = stack.front();
                stack.pop_front();
                if (update_rank(n)) {
                    nout = nei_out(n);
                    stack.insert(stack.end(), nout.begin(), nout.end());
                }
            }
        }

    void
        update_all_ranks()
        {
            std::deque<node_index_type> stack;
            std::vector<bool> visited(rank.size(), false);
            for (node_index_type n: active_nodes()) {
                if (nei_in(n).size() == 0) {
                    rank[n] = 0;
                    auto out = nei_out(n);
                    stack.insert(stack.end(), out.begin(), out.end());
                    visited[n] = true;
                }
            }
            while (stack.size()) {
                node_index_type n = stack.front();
                stack.pop_front();
                if (!visited[n]) {
                    update_rank(n);
                    auto out = nei_out(n);
                    stack.insert(stack.end(), out.begin(), out.end());
                    visited[n] = true;
                }
            }
        }

    void
        check_neighbours()
        {
            for (node_index_type i : active_nodes()) {
                for (node_index_type n: nei_in(i)) {
                    if (represented_by[n] != n) {
                        MSG_DEBUG("error neighbour_in[" << i << "] " << n << " is represented by " << represented_by[n]);
                    }
                    auto tmp = nei_out(n);
                    if (std::find(tmp.begin(), tmp.end(), i) == tmp.end()) {
                        MSG_DEBUG("error neighbour_in[" << i << "] " << n << " doesn't have " << i << " as a neighbour_out " << tmp);
                    }