graphnode.h 59.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#ifndef _SPEL_BAYES_GRAPH_NODE_H_
#define _SPEL_BAYES_GRAPH_NODE_H_

#include <memory>
#include <vector>
#include <algorithm>
#include <iostream>
#include <deque>
#include <list>
#include <map>
11
#include <cstdio>
12

13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
/*#include "factor.h"*/
#include "../pedigree.h"

template <typename V>
std::vector<V>
operator + (const std::vector<V>& u, const std::vector<V>& v)
{
    std::vector<V> ret(u.size() + v.size());
    auto it = std::set_union(u.begin(), u.end(), v.begin(), v.end(), ret.begin());
    ret.resize(it - ret.begin());
    return ret;
}


template <typename V>
std::vector<V>
operator - (const std::vector<V>& u, const std::vector<V>& v)
{
    std::vector<V> ret(u.size() + v.size());
    auto it = std::set_difference(u.begin(), u.end(), v.begin(), v.end(), ret.begin());
    ret.resize(it - ret.begin());
    return ret;
}


template <typename V>
std::vector<V>
operator % (const std::vector<V>& u, const std::vector<V>& v)
{
    std::vector<V> ret(std::min(u.size(), v.size()));
    auto it = std::set_intersection(u.begin(), u.end(), v.begin(), v.end(), ret.begin());
    ret.resize(it - ret.begin());
    return ret;
}




typedef int variable_index_type;
typedef size_t node_index_type;
typedef std::vector<node_index_type> node_vec;
typedef std::vector<variable_index_type> var_vec;
struct graph_type;
struct edge_type {
    const graph_type* graph;
    node_index_type first, second;

    edge_type(const graph_type* g, node_index_type f, node_index_type s) : graph(g), first(f), second(s) {}

    bool
        operator < (const edge_type& other) const
        {
            return graph < other.graph
                || (graph == other.graph
                        && (first < other.first
                            || (first == other.first && second < other.second)));
        }
};


#ifndef MESSAGE

static inline std::string __fetch_string(const std::ostream& os)
{
    return dynamic_cast<const std::stringstream*>(&os)->str();
}

#define MESSAGE(_expr_) __fetch_string(std::stringstream() << _expr_)

#endif
83
84
85
86
87
88


struct colour_proxy_impl;
typedef std::shared_ptr<colour_proxy_impl> colour_proxy;
struct colour_proxy_impl {
    colour_proxy proxy;
89
    colour_proxy cache;
90
91
92
93
94

    friend
        colour_proxy
        get_colour_impl(colour_proxy col)
        {
95
96
            while (col->cache->proxy) { col->cache = col->cache->proxy; }
            return col->cache;
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
        }

    friend
        colour_proxy
        assign_colour_impl(colour_proxy old_col, colour_proxy new_col)
        {
            get_colour_impl(old_col)->proxy = get_colour_impl(new_col);
            return old_col;
        }

    friend
        bool
        colour_equal(colour_proxy c1, colour_proxy c2)
        {
            return get_colour_impl(c1) == get_colour_impl(c2);
        }
};

inline
colour_proxy
117
create_colour() { auto ret = std::make_shared<colour_proxy_impl>(); ret->cache = ret; return ret; }
118
119


120
template <typename V>
121
void
122
sort_and_unique(std::vector<V>& v)
123
124
125
126
127
128
{
    std::sort(v.begin(), v.end());
    v.erase(std::unique(v.begin(), v.end()), v.end());
}


129
enum node_type { Factor, Interface, Aggregate };
130

131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330

template <typename V>
std::ostream&
operator << (std::ostream& os, const std::list<V>& v)
{
    auto i = v.begin(), j = v.end();
    if (i != j) {
        os << (*i);
        for (++i; i != j; ++i) {
            os << " -- " << (*i);
        }
    }
    return os;
}


template <typename K, typename V>
std::ostream&
operator << (std::ostream& os, const std::map<K, V>& m)
{
    os << "{ ";
    for (const auto& kv: m) {
        os << kv.first << ": " << kv.second << ", ";
    }
    return os << '}';
}




typedef std::vector<genotype_comb_type> message_type;


inline
message_type&
operator *= (message_type& accum, const message_type& msg)
{
    message_type tmp;

    std::map<variable_index_type, colour_proxy> colours;

    std::vector<var_vec> accum_varsets, msg_varsets;

    for (size_t ai = 0; ai < accum.size(); ++ai) {
        accum_varsets.push_back(get_parents(accum[ai]));
        for (variable_index_type v: accum_varsets.back()) {
            auto& ptr = colours[v];
            if (!ptr) {
                ptr = create_colour();
            }
        }
    }
    for (size_t mi = 0; mi < msg.size(); ++mi) {
        msg_varsets.push_back(get_parents(msg[mi]));
        for (variable_index_type v: msg_varsets.back()) {
            auto& ptr = colours[v];
            if (!ptr) {
                ptr = create_colour();
            }
        }
    }

    for (size_t i = 0; i < accum.size(); ++i) {
        auto mcol = colours[accum_varsets[i].front()];
        for (variable_index_type v: accum_varsets[i]) {
            auto& vcol = colours[v];
            if (vcol && !colour_equal(vcol, mcol)) {
                assign_colour_impl(vcol, mcol);
            }
        }
    }

    for (size_t i = 0; i < msg.size(); ++i) {
        auto mcol = colours[msg_varsets[i].front()];
        for (variable_index_type v: msg_varsets[i]) {
            auto& vcol = colours[v];
            if (vcol && !colour_equal(vcol, mcol)) {
                assign_colour_impl(vcol, mcol);
            }
        }
    }

    std::vector<colour_proxy> uniq;
    for (const auto& kv: colours) {
        uniq.push_back(get_colour_impl(kv.second));
    }
    sort_and_unique(uniq);

    std::vector<size_t> accum_dest, msg_dest;
    for (size_t i = 0; i < accum.size(); ++i) {
        accum_dest.push_back(std::find(uniq.begin(), uniq.end(), get_colour_impl(colours[accum_varsets[i].front()])) - uniq.begin());
    }
    for (size_t i = 0; i < msg.size(); ++i) {
        msg_dest.push_back(std::find(uniq.begin(), uniq.end(), get_colour_impl(colours[msg_varsets[i].front()])) - uniq.begin());
    }

    tmp.resize(uniq.size());

    for (size_t i = 0; i < accum.size(); ++i) {
        auto& table = tmp[accum_dest[i]];
        if (table.size()) {
            table = kronecker(table, accum[i]);
        } else {
            table = accum[i];
        }
    }
    for (size_t i = 0; i < msg.size(); ++i) {
        auto& table = tmp[msg_dest[i]];
        if (table.size()) {
            table = kronecker(table, msg[i]);
        } else {
            table = msg[i];
        }
    }

    tmp.swap(accum);

    return accum;
}

inline
std::shared_ptr<message_type>
operator *= (std::shared_ptr<message_type> accum, const std::shared_ptr<message_type>& msg)
{
    *accum *= *msg;
    return accum;
}


inline
message_type
operator % (const message_type& msg, const var_vec& variables)
{
    message_type tmp;
    tmp.reserve(msg.size());
    for (const auto& table: msg) {
        auto varset = get_parents(table);
        if (varset == variables) {
            tmp.push_back(table);
        } else {
            auto proj = varset % variables;
            if (proj.size()) {
                auto norm = varset - proj;
                tmp.push_back(project(table, proj, norm));
            }
        }
    }
    return tmp;
}


inline
std::shared_ptr<message_type>
operator % (std::shared_ptr<message_type> msg, const var_vec& variables)
{
    return std::make_shared<message_type>(*msg % variables);
}


inline
std::shared_ptr<message_type>
operator %= (std::shared_ptr<message_type> msg, const var_vec& variables)
{
    auto tmp = (*msg) % variables;
    msg->swap(tmp);
    return msg;
}


inline
genotype_comb_type&
operator *= (genotype_comb_type& table, const message_type& msg)
{
    if (table.size()) {
        for (const auto& mt: msg) {
            table = kronecker(table, mt);
        }
    } else {
        auto i = msg.begin(), j = msg.end();
        if (i != j) {
            table = *i++;
        }
        for (; i != j; ++i) {
            table = kronecker(table, *i);
        }
    }
    return table;
}












331
struct graph_type {
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419

    struct query_op_atom_type {
        const graph_type& g;
        var_vec variables;
        node_index_type node;
        size_t n_operands;

        query_op_atom_type(const graph_type& _g, const var_vec& vv, node_index_type n, size_t no)
            : g(_g), variables(vv), node(n), n_operands(no)
        {}

        friend
            std::ostream&
            operator << (std::ostream& os, const query_op_atom_type& qoa)
            {
                return os << "<op out=" << qoa.variables << " inner=" << (qoa.g.is_interface(qoa.node) ? 'I' : 'F') << qoa.g.variables_of(qoa.node) << " n_opd=" << qoa.n_operands << '>';
            }

        void
            operator () (std::vector<std::shared_ptr<message_type>>& stack, std::vector<const var_vec*>& var_stack) const
            {
                if (g.is_interface(node)) {
                    stack.emplace_back(std::make_shared<message_type>());
                } else {
                    auto ret = std::make_shared<message_type>();
                    for (size_t i = 0; i < n_operands; ++i) {
                        auto op = stack.back();
                        stack.pop_back();
                        var_stack.pop_back();
                        *ret *= *op;
                    }
                    ret %= variables;
                    stack.emplace_back(ret);
                }
                var_stack.push_back(&variables);
            }
    };


    typedef std::vector<query_op_atom_type> query_operation_type;

    enum ComputeStateOperation { PushFactor, PushMessage, Accumulate, Project, Store };
    struct compute_state_operation_type {
        const graph_type& g;
        ComputeStateOperation op_type;
        size_t n;
        edge_type e;
        var_vec v;
        query_operation_type op;

        compute_state_operation_type(const graph_type& _g, ComputeStateOperation ot, size_t _n, const edge_type& _e, var_vec&& _v)
            : g(_g), op_type(ot), n(_n), e(_e), v(std::move(_v))
        {}

        void
            operator () (std::map<edge_type, std::shared_ptr<message_type>>& messages, std::vector<std::shared_ptr<message_type>>& stack) const
            {
                switch (op_type) {
                    case PushFactor:
                        if (g.is_aggregate(n)) {
                            g.subgraphs[n]->compute_state(stack.end() - n, stack.end());
                            stack.push_back(g.subgraphs[n]->extract(op));
                        } else {
                            stack.push_back(g.tables[n]);
                        }
                        break;
                    case PushMessage:
                        stack.push_back(messages[e]);
                        break;
                    case Accumulate:
                        for (size_t i = 0; i < n; ++i) {
                            auto m2 = stack.back();
                            stack.pop_back();
                            stack.back() *= m2;
                        }
                        break;
                    case Project:
                        stack.back() %= v;
                        break;
                    case Store:
                        messages[e] = stack.back();
                        stack.pop_back();
                        break;
                    default:;
                };
            }
    };

420
421
422
423
424
425
426
427
428
429

    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;
430

431
432
    std::map<variable_index_type, node_index_type> interface_to_node;
    std::map<node_index_type, variable_index_type> node_to_interface;
433

434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
    std::vector<std::shared_ptr<graph_type>> subgraphs;

    std::vector<std::shared_ptr<message_type>> tables;
    std::vector<message_type> state;
    std::map<variable_index_type, std::vector<bn_label_type>> 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 = 1;

    graph_type()
        : rank(), represented_by(), type(), colour(), neighbours_in(), neighbours_out(), inner_nodes(), rules(), variables(), interface_to_node(), node_to_interface(), subgraphs(),  parent(nullptr), index_in_parent(0),aggregate_cycles(true), generate_interfaces(true), n_alleles(1)
    {}

    graph_type(size_t n_al)
        : rank(), represented_by(), type(), colour(), neighbours_in(), neighbours_out(), inner_nodes(), rules(), variables(), interface_to_node(), node_to_interface(), subgraphs(), parent(nullptr), index_in_parent(0), aggregate_cycles(true), generate_interfaces(true), n_alleles(n_al)
    {}

458

459
    bool is_aggregate(node_index_type node) const { return inner_nodes[node].size() > 1; }
460

461
462
463
464
465
466
467
468
469
470
471
472
473
    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;
        }
474

475
476
477
478
479
480
481
482
483
484
485
486
487
    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;
        }
488

489
490
491
492
493
    bool
        is_computable(node_index_type node) const
        {
            return !is_interface(node);
        }
494
495
496

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

497
    node_vec
498
499
        active_nodes() const
        {
500
            node_vec ret;
501
            ret.reserve(represented_by.size());
502
            for (node_index_type i = 0; i < represented_by.size(); ++i) {
503
504
505
506
507
508
509
                if (represented_by[i] == i) {
                    ret.push_back(i);
                }
            }
            return ret;
        }

510
511
    node_vec
        resolve_vector(const node_vec& vec) const
512
        {
513
            node_vec ret;
514
            ret.reserve(vec.size());
515
            for (node_index_type i: vec) {
516
517
518
519
520
521
                ret.push_back(resolve(i));
            }
            sort_and_unique(ret);
            return ret;
        }

522
523
    var_vec
        interface_nodes(var_vec inputs) const
524
        {
525
526
527
528
529
530
531
532
            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;
        }

533
    std::vector<edge_type>
534
535
        active_edges() const
        {
536
            std::vector<edge_type> ret;
537
538
            for (node_index_type n: active_nodes()) {
                for (node_index_type o: resolve_vector(neighbours_out[n])) {
539
                    ret.emplace_back(this, n, o);
540
541
542
543
544
545
546
547
                }
            }
            return ret;
        }


    void
        dump_node(node_index_type n)
548
549
        {
            std::cout
550
                << '[' << rank[n] << "] " << (is_interface(n) ? "INTERFACE " : (n == inner_nodes[n][0] ? "FACTOR " : "AGGREGATE ")) << n << std::endl
551
552
553
554
555
556
                << "  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
557
                << "  variable(s) " << variables_of(n) << std::endl
558
559
                ;
            if (inner_nodes[n].size() > 1) {
560
561
562
563
564
                for (node_index_type i: inner_nodes[n]) {
                    /*var_vec remaining, imported;*/
                    /*restrict_inputs(rules[i], inner_nodes[n], remaining, imported);*/
                    /*std::cout << "  * rule for " << i << ": " << remaining << " / " << imported << std::endl;*/
                    std::cout << "  * rule for " << i << ": " << rules[i] << std::endl;
565
566
567
568
569
570
                }
            }
            std::cout
                << std::endl;
        }

571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
#define DUMP_SZ(_x) << " * " #_x " " << _x.size() << std::endl

    void
        dump_sizes() const
        {
            std::cout
                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)
                ;
        }

589
590
591
592
    void
        dump()
        {
            std::cout << "ALL NODES" << std::endl;
593
594
            dump_sizes();
            for (node_index_type i = 0; i < rank.size(); ++i) {
595
596
597
598
599
600
601
602
                dump_node(i);
            }
        }

    void
        dump_active()
        {
            std::cout << "ACTIVE NODES" << std::endl;
603
604
            dump_sizes();
            for (node_index_type i: active_nodes()) {
605
606
607
608
                dump_node(i);
            }
        }

609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
    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 = resolve_vector(neighbours_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];
                }
            }
        }

    void
637
        finalize()
638
        {
639
640
641
642
643
            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); }
            };
            /*dump_active();*/
644
            auto edges = active_edges();
645
            std::map<emitter_and_interface_type, size_t, compare_eai> interface_map;
646
647
            for (const auto& e: edges) {
                if (is_interface(e.first) ^ !is_interface(e.second)) {
648
                    /* By construction, it's a factor->factor edge, never an interface->interface edge */
649
                    auto varset = variables_of(e.first) % variables_of(e.second);
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
                    /*std::cout << "edge between factor(graph)s " << e.first << " and " << e.second << " carries variables " << varset << std::endl;*/
                    node_index_type& i = interface_map[{e.first, varset}];
                    if (i == 0) {
                        if (varset.size() == 1) {
                            i = add_interface(node_vec{e.first}, varset.front());
                            colour[i] = get_colour_impl(colour[e.first]);
                            neighbours_out[i].push_back(e.second);
                        } 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{e.first}, node_vec{e.second}, var_vec{}, get_colour_impl(colour[e.first]), Aggregate, iv, -1);
                            for (node_index_type ni: iv) {
                                represented_by[ni] = i;
                            }
667
                        }
668
669
                    } else {
                        neighbours_out[i].push_back(e.second);
670
671
672
                    }
                    filter_out_and_replace_by(neighbours_out[e.first], node_vec{e.second}, i);
                    filter_out_and_replace_by(neighbours_in[e.second], node_vec{e.first}, i);
673
674
675
                    /*dump_node(e.first);*/
                    /*dump_node(i);*/
                    /*dump_node(e.second);*/
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
                }
            }

            if (size()) {
                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;
                        }
                    }
                }

                for (node_index_type n: A) {
                    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);
                            neighbours_out[n].push_back(i);
                        }
                    }
                }
700
701
702
703
704
705
706

                for (node_index_type n: A) {
                    if (is_aggregate(n) && !is_interface(n)) {
                        std::cout << "subgraphs.size() = " << subgraphs.size() << " n = " << n << std::endl;
                        subgraphs[n] = subgraph(n);
                    }
                }
707
708
709
710
711
712
713
            }
        }

    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)
714
        {
715
            node_index_type ret = rank.size();
716
            /*std::cout << "adding node " << ret << std::endl;*/
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
            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);
            }
            /*std::cout << "rank=" << rank.back() << std::endl;*/
            type.push_back(t);
            variables.push_back(var);
733
            represented_by.push_back(ret);
734
735
            subgraphs.emplace_back();
            tables.emplace_back();
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
            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()) {
                /*std::cout << "resolve_interface(" << var << ") => new interface" << std::endl;*/
                return add_interface(node_vec{}, var);
            }
            /*std::cout << "resolve_interface(" << var << ") => resolve(" << it->second << ") = " << resolve(it->second) << std::endl;*/
            return resolve(it->second);
        }

    node_index_type
        add_interface(const node_vec& producer, variable_index_type var)
        {
            node_index_type ret = add_node(producer, node_vec{}, var_vec{}, create_colour(), Interface, node_vec{}, var);
            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)
        {
806
            /*std::cout << "add_factor(" << var << ')' << std::endl;*/
807
808
            node_index_type ret = add_factor(var_vec{}, create_colour(), var);
            /*dump();*/
809
810
811
            return ret;
        }

812
813
    node_index_type
        resolve(node_index_type n) const
814
815
816
817
818
        {
            while (n != represented_by[n]) { n = represented_by[n]; }
            return n;
        }

819
820
    node_index_type
        add_factor(variable_index_type v1, variable_index_type var)
821
        {
822
            /*std::cout << "add_factor(" << v1 << ", " << var << ')' << std::endl;*/
823
824
825
826
827
828
829
830
831
832
833
834
            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);
835
            }
836
837
838
            node_index_type ret = add_factor(var_vec{v1}, colour[p1r], var);
            /*compute_ranks();*/
            /*dump();*/
839
840
841
            return ret;
        }

842
843
    node_index_type
        add_factor(variable_index_type v1,variable_index_type v2, variable_index_type var)
844
        {
845
            /*std::cout << "add_factor(" << v1 << ", " << v2 << ", " << var << ')' << std::endl;*/
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
            node_index_type p1r, p2r;
            if (!generate_interfaces) {
                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 {
                p1r = resolve_interface(v1);
                p2r = resolve_interface(v2);
            }

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

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

            bool cycle = false;
            if (p1r != p2r) {
                if (colour_equal(colour[p1r], colour[p2r])) {
                    /* cycle! */
881
                    /*std::cout << "cycle!" << std::endl;*/
882
883
                    cycle = true;
                } else {
884
                    /*std::cout << "no cycle." << std::endl;*/
885
886
887
888
889
                    assign_colour_impl(colour[p1r], colour[p2r]);
                }
            }

            if (cycle && aggregate_cycles) {
890
                /*std::cout << "search a path between " << p1r << " and " << p2r << std::endl;*/
891
892
                auto path = find_path_between_parents(p1r, p2r, ret);
                /*dump_active();*/
893
                /*std::cout << "Found path: "; for (size_t n: path) { std::cout << ' ' << n; } std::cout << std::endl;*/
894
                parents_of_max_rank aggr_first;
895
                /*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]; })];*/
896
                while (path.size() > 2 && ((aggr_first = find_parents_of_max_rank(path)), !aggregate(*aggr_first.p1, *aggr_first.p2, ret, path, aggr_first.p1))) {
897
898
                    path.erase(aggr_first.p1);
                    path.erase(aggr_first.p2);
899
900
901
902
                    path.erase(aggr_first.child);
                }
            }

903
904
905
906
            /*compute_ranks();*/
            /*dump();*/
            return ret;
        }
907

908
909
910
911
912
913
914
915
    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.emplace_back(letter, letter, al, al);
916
            }
917
        }
918

919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
    void
        compute_labels(genotype_comb_type& table, variable_index_type spawnling)
        {
            auto& domain = domains[spawnling];
            domain.clear();
            domain.reserve(table.m_combination.size());
            for (auto& state: table) {
                auto& keys = state.keys.keys;
                bn_label_type G = keys.front().state;
                const auto& p1 = keys[G.first_allele];
                const auto& p2 = keys[G.second_allele];
                bool f1 = G.first == GAMETE_L;
                bool f2 = G.second == GAMETE_L;
                /*std::cout << "keys " << keys << " G " << G << " p1 " << p1 << " p2 " << p2 << " f1 " << f1 << " f2 " << f2 << std::endl;*/
                domain.emplace_back(
                        f1 ? p1.state.first : p1.state.second,
                        f2 ? p2.state.first : p2.state.second,
                        f1 ? p1.state.first_allele : p1.state.second_allele,
                        f2 ? p2.state.first_allele : p2.state.second_allele
                        );
                keys.emplace_back(spawnling, domain.back());
940
            }
941
942
            sort_and_unique(domain);
        }
943

944
945
946
947
948
949
950
    std::shared_ptr<message_type>
        compute_factor_table(variable_index_type spawnling, const var_vec& parents, bool dh)
        {
            auto op = build_query_operation(parents);
            genotype_comb_type jpar;
            for (const auto& expr: op) {
                jpar *= *extract(expr);
951
            }
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
            std::cout << "Computed parent probability table " << jpar << std::endl;

            std::vector<bn_label_type> label_g;
            if (parents.size() == 2) {
                label_g = {
                    {GAMETE_L, GAMETE_L, 1, 2},
                    {GAMETE_L, GAMETE_R, 1, 2},
                    {GAMETE_R, GAMETE_L, 1, 2},
                    {GAMETE_R, GAMETE_R, 1, 2}
                };
            } else if (dh) {
                label_g = {
                    {GAMETE_L, GAMETE_L, 1, 1},
                    {GAMETE_R, GAMETE_L, 1, 1},
                };
            } else {
                label_g = {
                    {GAMETE_L, GAMETE_L, 1, 1},
                    {GAMETE_L, GAMETE_R, 1, 1},
                    {GAMETE_R, GAMETE_L, 1, 1},
                    {GAMETE_R, GAMETE_R, 1, 1}
                };
974
975
            }

976
977
            genotype_comb_type
                G = state_to_combination(-1, label_g) * (1. / label_g.size());
978

979
980
981
982
983
984
985
986
987
            jpar = kronecker(jpar, G);
            std::cout << "Computed parent probability table with gametes " << jpar << std::endl;
            compute_labels(jpar, spawnling);
            std::cout << "Computed factor table with gametes " << jpar << std::endl;
            std::cout << "Computed domain for #" << spawnling << ": " << domains[spawnling] << std::endl;
            auto ret = std::make_shared<message_type>();
            ret->emplace_back(fold(sum_over(jpar, {-1})));
            return ret;
        }
988

989
990
991
992
993
994
995
996
997
998
    node_index_type
        add_cross(variable_index_type p1, variable_index_type p2, variable_index_type id)
        {
            auto factor = compute_factor_table(id, var_vec{p1, p2}, false);
            node_index_type ret = add_factor(p1, p2, id);
            tables[ret] = factor;
            std::cout << "Computed factor for #" << ret << ": " << (*tables[ret]) << std::endl;
            /*operations.push_back(op);*/
            return ret;
        }
999

1000
    node_index_type