preprocess.h

Go to the documentation of this file.

// Copyright 2024 Markus Anders
// This file is part of dejavu 2.0.
// See LICENSE for extended copyright information.
 
#ifndef SASSY_H
#define SASSY_H
 
#define SASSY_VERSION_MAJOR 2
#define SASSY_VERSION_MINOR 0
 
#include "graph.h"
#include "ds.h"
#include "refinement.h"
#include "components.h"
#include <vector>
#include <iomanip>
#include <ctime>
 
namespace dejavu {
    using dejavu::ds::coloring;
 
    class preprocessor;
    static preprocessor* save_preprocessor;
 
    class preprocessor {
    private:
        enum preop { deg01, deg2ue, deg2ma, qcedgeflip, densify2, twins };
 
        //inline static preprocessor* save_preprocessor;
        dejavu_hook* saved_hook = nullptr;
 
        dejavu::timed_print* print = nullptr;
 
        coloring c;
        dejavu::worklist automorphism;
        dejavu::worklist automorphism_supp;
 
        dejavu::worklist aux_automorphism;
        dejavu::worklist aux_automorphism_supp;
 
        bool layers_melded = false;
 
        bool skipped_preprocessing = false;
 
        dejavu::markset del;
        dejavu::markset del_e;
 
        std::vector<std::vector<int>> translation_layers;
        std::vector<std::vector<int>> backward_translation_layers;
 
        std::vector<int>              backward_translation;
        std::vector<std::vector<int>> recovery_strings;
 
        std::vector<std::pair<int, int>> baked_recovery_string_pt;
        std::vector<int>                 baked_recovery_string;
 
        std::vector<std::vector<int>> add_edge_buff;
        dejavu::worklist worklist_deg0;
        dejavu::worklist worklist_deg1;
        dejavu::markset add_edge_buff_act;
        dejavu::ir::refinement* R1 = nullptr;
 
        dejavu::markset touched_color_cache;
 
        std::vector<int> g_old_v;
        std::vector<int> g_old_e;
        dejavu::worklist edge_scratch;
 
        std::vector<int> translate_layer_fwd;
        std::vector<int> translate_layer_bwd;
 
        int edge_attached_information = 0;
        std::vector<int>                                  recovery_edge_attached;
        std::vector<std::vector<std::tuple<int,int,int>>> recovery_edge_adjacent;
 
        dejavu::worklist _automorphism;
        dejavu::worklist _automorphism_supp;
        std::vector<int> save_colmap;
 
        dejavu::worklist before_move;
 
        dejavu::ir::graph_decomposer* decomposer = nullptr;
        int current_component = 0;
 
        bool ir_quotient_component_init = false;
 
    public:
        dejavu::big_number grp_sz; 
        int domain_size   = 0;      
        bool h_deact_deg1 = false;  
        bool h_deact_deg2 = false;  
        preprocessor() = default;
        explicit preprocessor(dejavu::ir::refinement* R) {
            R1 = R;
        }
 
        explicit preprocessor(dejavu::timed_print* printer) : print(printer) {};
 
 
        // for a vertex v of reduced graph, return corresponding vertex of the original graph
        int translate_back(int v) {
            const int layers = static_cast<int>(translation_layers.size());
            for (int l = layers - 1; l >= 0; --l) {
                v = backward_translation_layers[l][v];
            }
            return v;
        }
 
        // keeps track of the group size generated by automorphisms found so far
        void multiply_to_group_size(int n) {
            grp_sz.multiply(n);
        }
 
        void inject_decomposer(dejavu::ir::graph_decomposer* new_decomposer, int component) {
            decomposer = new_decomposer;
            current_component = component;
        }
 
    private:
        // combine translation layers
        void meld_translation_layers() {
            if (layers_melded)
                return;
            const int layers = static_cast<int>(translation_layers.size());
            if(layers == 0) {
                backward_translation_layers.emplace_back();
                backward_translation_layers[0].reserve(domain_size);
                for(int i = 0; i < domain_size; ++i) {
                    backward_translation_layers[0].push_back(i);
                }
                layers_melded = true;
                return;
            }
 
 
            const int reduced_size = static_cast<int>(backward_translation_layers[layers - 1].size());
            backward_translation.reserve(reduced_size);
            for (int i = 0; i < reduced_size; ++i) backward_translation.push_back(i);
            for (int i = 0; i < reduced_size; ++i) {
                int next_v = i;
                for (int l = layers - 1; l >= 0; --l) { // TODO inefficient
                    next_v = backward_translation_layers[l][next_v];
                }
                backward_translation[i] = next_v;
            }
            layers_melded = true;
 
            baked_recovery_string.reserve(domain_size);
            baked_recovery_string_pt.reserve(domain_size);
 
            for (int i = 0; i < reduced_size; ++i) {
                const int orig_i = backward_translation[i];
                const int pt_start = static_cast<int>(baked_recovery_string.size());
                for(auto v : recovery_strings[orig_i]) baked_recovery_string.push_back(v);
                const int pt_end = static_cast<int>(baked_recovery_string.size());
                baked_recovery_string_pt.emplace_back(pt_start, pt_end);
            }
            recovery_strings.clear();
        }
 
        int twins_partition_refinement(dejavu::sgraph *g, int *colmap, dejavu_hook* hook, bool false_twins) {
            g->initialize_coloring(&c, colmap);
            copy_coloring_to_colmap(&c, colmap);
            del.reset();
 
            dejavu::ds::worklist neighbour_colors(g->v_size);
            dejavu::ds::markset  is_neighbour(g->v_size);
            workset_t<int>       neighbours(g->v_size);
            workspace            scratch(2*g->v_size);
 
            int removed_twins = 0;
 
            for(int refine_with_v = 0; refine_with_v < g->v_size; ++refine_with_v) {
                const int start_pt = g->v[refine_with_v];
                const int end_pt   = start_pt + g->d[refine_with_v];
 
                neighbour_colors.reset();
                neighbours.reset();
 
                for(int j = start_pt; j < end_pt; ++j) {
                    const int neighbour     = g->e[j];
                    const int neighbour_col = c.vertex_to_col[neighbour];
 
                    // singletons can not be split, so let's just ignore them
                    if (c.ptn[neighbour_col] == 0) continue;
 
                    if(neighbours.get(neighbour_col) == -1) {
                        neighbour_colors.push_back(neighbour_col);
                        // neighbours acts as the write position for degree 1 vertices of col in the scratchpad
                        neighbours.set(neighbour_col, 0);
                    }
 
                    // write vertex to scratchpad for later use
                    dej_assert(2*neighbour_col + neighbours.get(neighbour_col) < 2*g->v_size);
                    scratch[2*neighbour_col + neighbours.get(neighbour_col)] = neighbour;
                    neighbours.inc_nr(neighbour_col); // we reset neighbours later, use old_color_classes for reset
                }
 
                if(false_twins) {
                    // do same as above, but for myself
                    const int neighbour     = refine_with_v;
                    const int neighbour_col = c.vertex_to_col[neighbour];
 
                    // singletons can not be split, so let's just ignore them
                    if (c.ptn[neighbour_col] > 0) {
                        if (neighbours.get(neighbour_col) == -1) {
                            neighbour_colors.push_back(neighbour_col);
                            // neighbours acts as the write position for degree 1 vertices of col in the scratchpad
                            neighbours.set(neighbour_col, 0);
                        }
                        // write vertex to scratchpad for later use
                        dej_assert(2 * neighbour_col + neighbours.get(neighbour_col) < 2 * g->v_size);
                        scratch[2 * neighbour_col + neighbours.get(neighbour_col)] = neighbour;
                        neighbours.inc_nr(neighbour_col);// we reset neighbours later, use old_color_classes for reset
                    }
                }
 
                for(int j = 0; j < neighbour_colors.cur_pos; ++j) {
                    const int deg0_col    = neighbour_colors[j];
                    const int deg1_col_sz = neighbours.get(deg0_col); //  - deg0_col
                    const int deg0_col_sz = (c.ptn[deg0_col] + 1) - deg1_col_sz;
                    const int deg1_col = deg0_col + deg0_col_sz;
 
 
                    dej_assert(c.vertex_to_col[c.lab[deg0_col]] == deg0_col);
 
                    // no split? done...
                    if (deg0_col == deg1_col) {
                        neighbours.set(deg0_col, -1);
                        dej_assert(c.vertex_to_col[c.lab[deg0_col]] == deg0_col);
                        continue;
                    }
 
                    // there is a split
                    ++c.cells;
 
                    // set ptn
                    c.ptn[deg0_col] = deg0_col_sz - 1;
                    c.ptn[deg1_col] = deg1_col_sz - 1;
                    c.ptn[deg1_col - 1] = 0;
 
                    int deg1_write_pos = deg1_col;
                    int deg1_read_pos = 2*deg0_col + neighbours.get(deg0_col) - 1;
 
                    // rearrange vertices of deg1 to the back of deg0 color
                    dej_assert(deg1_read_pos >= 2*deg0_col);
 
                    while (deg1_read_pos >= 2*deg0_col) {
                        dej_assert(deg1_read_pos < 2*g->v_size);
                        const int v = scratch[deg1_read_pos];
                        const int vertex_at_pos = c.lab[deg1_write_pos];
                        const int lab_pos = c.vertex_to_lab[v];
 
                        c.lab[deg1_write_pos] = v;
                        c.vertex_to_lab[v] = deg1_write_pos;
                        c.vertex_to_col[v] = deg1_col;
 
                        c.lab[lab_pos] = vertex_at_pos;
                        c.vertex_to_lab[vertex_at_pos] = lab_pos;
 
                        deg1_write_pos++;
                        deg1_read_pos--;
                    }
 
                    dej_assert(deg1_write_pos == deg1_col + deg1_col_sz);
                    dej_assert(c.vertex_to_col[c.lab[deg0_col]] == deg0_col);
                    dej_assert(c.vertex_to_col[c.lab[deg1_col]] == deg1_col);
 
                    // reset neighbours count to -1
                    neighbours.set(deg0_col, -1);
                }
 
                if(c.cells == g->v_size) return removed_twins;
            }
 
            for(int i = 0; i < g->v_size;) {
                const int twin_class    = i;
                const int twin_class_sz = c.ptn[twin_class] + 1;
                i += twin_class_sz;
 
                if(twin_class_sz > 1) {
                    const int remaining_vertex = c.lab[twin_class];
                    const int orig_remaining_vertex = translate_back(remaining_vertex);
                    int group_size_calculation = 2;
 
                    // remove all but one twin of this class and output transposition automorphisms
                    for (int j = twin_class + 1; j < twin_class + twin_class_sz; ++j) {
                        const int vertex = c.lab[j];
 
                        // multiply group size
                        multiply_to_group_size(group_size_calculation);
                        ++group_size_calculation;
 
                        // transposition of remaining_vertex and vertex
                        _automorphism[vertex]           = remaining_vertex;
                        _automorphism[remaining_vertex] = vertex;
 
                        dej_assert(colmap[vertex] == colmap[remaining_vertex]);
 
                        _automorphism_supp.push_back(vertex);
                        _automorphism_supp.push_back(remaining_vertex);
 
                        pre_hook(g->v_size, _automorphism.get_array(), _automorphism_supp.cur_pos,
                                 _automorphism_supp.get_array(), hook);
 
                        reset_automorphism(_automorphism.get_array(), _automorphism_supp.cur_pos,
                                           _automorphism_supp.get_array());
                        _automorphism_supp.reset();
                    }
 
                    for (int j = twin_class + 1; j < twin_class + twin_class_sz; ++j) {
                        const int vertex = c.lab[j];
 
                        // add to recovery string of remaining vertex
                        const int orig_other = translate_back(vertex);
                        recovery_strings[orig_remaining_vertex].push_back(orig_other);
                        for (int k = 0; k < static_cast<int>(recovery_strings[orig_other].size()); ++k) {
                            recovery_strings[orig_remaining_vertex].push_back(recovery_strings[orig_other][k]);
                        }
 
                        // remove vertex later
                        del.set(vertex);
                        ++removed_twins;
                    }
 
                    // color the remaining vertex appropriately with the size of the twin class
                    colmap[remaining_vertex] = colmap[remaining_vertex] + (twin_class_sz-1);
                }
            }
 
            return removed_twins;
        }
 
        int remove_twins(dejavu::sgraph *g, int *colmap, dejavu_hook* hook) {
            //coloring col;
            g->initialize_coloring(&c, colmap);
 
            dejavu::ds::markset  test_twin(g->v_size);
            std::vector<int> add_to_string;
            dejavu::ds::worklist twin_counter(g->v_size);
 
            dejavu::ds::markset  potential_twin(g->v_size);
            dejavu::ds::worklist potential_twin_counter(g->v_size);
 
            dejavu::ds::markset  touched(g->v_size);
 
            std::vector<int> potential_twin_list;
 
            del.reset();
 
            int d_limit = std::max(static_cast<int>(sqrt(sqrt(g->v_size))), 8);
 
            int s_twins = 0;
            for(int i = 0; i < g->v_size; ++i) twin_counter[i] = 0;
                // iterate over color classes
            for(int i = 0; i < g->v_size;) {
                const int color = i;
                const int color_size = c.ptn[color] + 1;
                i += color_size;
 
                for(int j = 0; j < color_size; ++j) {
                    if(j == 2 && s_twins == 0) break;
                    const int vertex = c.lab[color + j];
                    if(g->d[vertex] > d_limit) break;
 
                    touched.set(vertex);
                    int twin_class_sz = 1;
 
                    bool limit_breached = false;
 
                    add_to_string.clear();
                    if(del.get(vertex)) continue;
                    test_twin.reset();
                    potential_twin.reset();
                    potential_twin_list.clear();
                    for(int k = 0; k < g->d[vertex]; ++k) {
                        const int neighbour = g->e[g->v[vertex] + k];
                        test_twin.set(neighbour);
                        if(g->d[neighbour] > d_limit) {
                            limit_breached = true;
                            break;
                        }
                        for(int l = 0; l < g->d[neighbour]; ++l) {
                            const int neighbour_neighbour = g->e[g->v[neighbour] + l];
                            if(touched.get(neighbour_neighbour)) continue;
                            if(potential_twin.get(neighbour_neighbour)) {
                                potential_twin_counter[neighbour_neighbour] += 1;
                            } else {
                                potential_twin.set(neighbour_neighbour);
                                potential_twin_counter[neighbour_neighbour] = 1;
                                potential_twin_list.push_back(neighbour_neighbour);
                            }
                        }
                    }
                    if(limit_breached) break;
 
                    //for(int jj = j+1; jj < color_size; ++jj) {
                    for(int other_vertex : potential_twin_list) {
                        bool is_twin = true;
                        //const int other_vertex = col.lab[color + jj];
                        if(potential_twin_counter[other_vertex] < g->d[vertex] - 1) continue;
                        if(vertex == other_vertex || del.get(other_vertex) ||
                                c.vertex_to_col[other_vertex] != c.vertex_to_col[vertex]) continue;
                        dej_assert(g->d[vertex] == g->d[other_vertex]);
                        for(int k = 0; k < g->d[other_vertex]; ++k) {
                            const int neighbour = g->e[g->v[other_vertex] + k];
                            if(neighbour == vertex) continue;
                            if(!test_twin.get(neighbour)) {
                                is_twin = false;
                                break;
                            }
                        }
 
                        if(is_twin) {
                            ++s_twins;
                            ++twin_counter[vertex];
                            dej_assert(vertex != other_vertex);
                            ++twin_class_sz;
                            multiply_to_group_size(twin_class_sz);
                            del.set(other_vertex);
                            add_to_string.push_back(other_vertex);
 
                            _automorphism[vertex]       = other_vertex;
                            _automorphism[other_vertex] = vertex;
                            _automorphism_supp.push_back(vertex);
                            _automorphism_supp.push_back(other_vertex);
 
                            pre_hook(g->v_size, _automorphism.get_array(), _automorphism_supp.cur_pos,
                                     _automorphism_supp.get_array(), hook);
 
                            reset_automorphism(_automorphism.get_array(), _automorphism_supp.cur_pos,
                                               _automorphism_supp.get_array());
                            _automorphism_supp.reset();
                        }
                    }
 
                    dej_assert(!del.get(vertex));
                    dej_assert(twin_class_sz -1 == static_cast<int>(add_to_string.size()));
                    dej_assert(twin_counter[vertex] == static_cast<int>(add_to_string.size()));
                    const int orig_vertex = translate_back(vertex);
                    for(auto& other_vertex : add_to_string) {
                        dej_assert(del.get(other_vertex));
                        const int orig_other  = translate_back(other_vertex);
                        dej_assert(static_cast<int>(recovery_strings[orig_vertex].size()) ==
                                   static_cast<int>(recovery_strings[orig_other].size()));
                        recovery_strings[orig_vertex].push_back(orig_other);
                        for(int k = 0; k < static_cast<int>(recovery_strings[orig_other].size()); ++k) {
                            recovery_strings[orig_vertex].push_back(recovery_strings[orig_other][k]);
                        }
                    }
                }
 
                add_to_string.clear();
                for(int j = 0; j < color_size; ++j) {
                    const int vertex = c.lab[color + j];
                    if(!del.get(vertex)) {
                        dej_assert(colmap[vertex] + twin_counter[vertex] < color + color_size);
                        colmap[vertex] = color + twin_counter[vertex];
                        add_to_string.push_back(vertex);
                    }
                }
            }
 
            return s_twins;
        }
 
        // reset internal automorphism structure to the identity
        static void reset_automorphism(int *rautomorphism, int nsupp, const int *supp) {
            for (int i = 0; i < nsupp; ++i) {
                rautomorphism[supp[i]] = supp[i];
            }
        };
 
        // are vert1 and vert2 adjacent in g?
        static bool is_adjacent(dejavu::sgraph*g, int vert1, int vert2) {
            if(g->d[vert1] < g->d[vert2]) {
                for(int i = 0; i < g->d[vert1]; ++i) {
                    if(g->e[g->v[vert1] + i] == vert2)
                        return true;
                }
                return false;
            } else {
                for(int i = 0; i < g->d[vert2]; ++i) {
                    if(g->e[g->v[vert2] + i] == vert1)
                        return true;
                }
                return false;
            }
        }
 
        // check for degree 2 matchings between color classes and mark for deletion
        void red_deg2_path_size_1(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1 || h_deact_deg2)
                return;
 
            del.reset();
 
            //coloring test_col;
            g->initialize_coloring(&c, colmap);
 
            assure_ir_quotient_init(g);
 
            touched_color_cache.reset();
 
            worklist_deg0.reset();
            worklist_deg1.reset();
 
            for (int i = 0; i < g->v_size; ++i) {
                worklist_deg1.push_back(-1);
                worklist_deg0.push_back(0);
            }
 
            add_edge_buff_act.reset();
            //add_edge_buff.reserve(g->v_size);
            //for (int i = 0; i < g->v_size; ++i)
             //   add_edge_buff.emplace_back(); // could do this smarter... i know how many edges end up here
 
            for (int i = 0; i < g->v_size;) {
                const int test_v = c.lab[i];
                const int path_col_sz = c.ptn[i];
                if (g->d[test_v] == 2) {
                    const int n1 = g->e[g->v[test_v] + 0];
                    const int n2 = g->e[g->v[test_v] + 1];
                    if (g->d[n1] != 2 && g->d[n2] != 2) {
                        // relevant path of length 1
                        const int col_n1 = c.vertex_to_col[n1];
                        const int col_n2 = c.vertex_to_col[n2];
 
                        const int col_sz_n1 = c.ptn[c.vertex_to_col[n1]];
                        const int col_sz_n2 = c.ptn[c.vertex_to_col[n2]];
 
                        if (col_sz_n1 != col_sz_n2 || col_sz_n1 != path_col_sz || col_n1 == col_n2) {
                            i += c.ptn[i] + 1;
                            continue;
                        }
 
                        bool already_matched_n1_n2 = false;
 
                        const int already_match_pt1 = g->v[c.lab[col_n1]];
                        const int already_match_pt2 = g->v[c.lab[col_n2]];
 
                        if (touched_color_cache.get(col_n1) && touched_color_cache.get(col_n2)) {
                            for (int j = 0; j < worklist_deg0[col_n1]; ++j) {
                                if (edge_scratch[already_match_pt1 + j] == col_n2)
                                    already_matched_n1_n2 = true;
                            }
                        }
 
                        if (touched_color_cache.get(col_n1) && touched_color_cache.get(col_n2) &&
                            already_matched_n1_n2 && worklist_deg0[col_n1] == 1 && worklist_deg0[col_n2] == 1) {
                            // const bool matching_same_cols1 =
                            //         test_col.vertex_to_col[worklist_deg1[n1]] == test_col.vertex_to_col[n2];
                            // const bool matching_same_cols2 =
                            //        test_col.vertex_to_col[worklist_deg1[n2]] == test_col.vertex_to_col[n1];
                            // const bool match_same_cols = matching_same_cols1 && matching_same_cols2;
 
                            bool check_if_match = true;
                            for (int f = 0; f < c.ptn[i] + 1; ++f) {
                                const int _test_v = c.lab[i + f];
                                const int _n1 = g->e[g->v[_test_v] + 0];
                                const int _n2 = g->e[g->v[_test_v] + 1];
                                check_if_match = check_if_match && (worklist_deg1[_n1] == _n2);
                                if (!check_if_match)
                                    break;
                            }
 
                            if (check_if_match) {
                                for (int f = 0; f < c.ptn[i] + 1; ++f) {
                                    const int _test_v = c.lab[i + f];
                                    del.set(_test_v);
                                }
                                for (int f = 0; f < c.ptn[i] + 1; ++f) {
                                    const int _test_v = c.lab[i + f];
                                    const int _n1 = g->e[g->v[_test_v] + 0];
                                    const int _n2 = g->e[g->v[_test_v] + 1];
                                    int can_n;
                                    if (c.vertex_to_col[_n1] < c.vertex_to_col[_n2])
                                        can_n = _n1;
                                    else
                                        can_n = _n2;
                                    const int orig_test_v = translate_back(_test_v);
                                    const int orig_n1 = translate_back(can_n);
                                    recovery_strings[orig_n1].push_back(orig_test_v);
                                    for (size_t s = 0; s < recovery_strings[orig_test_v].size(); ++s)
                                        recovery_strings[orig_n1].push_back(
                                                recovery_strings[orig_test_v][s]);
                                }
                            }
 
 
                            i += c.ptn[i] + 1;
                            continue;
                        }
 
                        if (touched_color_cache.get(col_n1) && touched_color_cache.get(col_n2) &&
                            already_matched_n1_n2) {
                            i += c.ptn[i] + 1;
                            continue;
                        }
 
                        const int col_endpoint2 = colmap[n2]; // colmap?
                        bool col_cycle = false;
                        for (int f = 0; f < g->d[n1]; ++f) {
                            const int col_other = colmap[g->e[g->v[n1] + f]];
                            if (col_other == col_endpoint2) {
                                col_cycle = true;
                                break;
                            }
                        }
 
                        if (col_cycle) {
                            i += c.ptn[i] + 1;
                            continue;
                        }
 
                        edge_scratch[already_match_pt1 + worklist_deg0[col_n1]] = col_n2;
                        // overwrites itself, need canonical vertex for color
                        edge_scratch[already_match_pt2 + worklist_deg0[col_n2]] = col_n1;
                        ++worklist_deg0[col_n1];
                        ++worklist_deg0[col_n2];
 
                        touched_color_cache.set(col_n1);
                        touched_color_cache.set(col_n2);
 
                        //const size_t debug_str_sz = recovery_strings[translate_back(test_v)].size();
 
                        for (int f = 0; f < c.ptn[i] + 1; ++f) {
                            const int _test_v = c.lab[i + f];
                            dej_assert(g->d[_test_v] == 2);
                            const int _n1 = g->e[g->v[_test_v] + 0];
                            const int _n2 = g->e[g->v[_test_v] + 1];
                            worklist_deg1[_n1] = _n2;
                            worklist_deg1[_n2] = _n1;
                            for (size_t t = 0; t < add_edge_buff[_n2].size(); ++t)
                                dej_assert(add_edge_buff[_n2][t] != _n1);
                            add_edge_buff[_n2].push_back(_n1);
                            add_edge_buff_act.set(_n2);
                            add_edge_buff[_n1].push_back(_n2);
                            add_edge_buff_act.set(_n1);
                            del.set(_test_v);
 
                            int can_n;
                            if (c.vertex_to_col[_n1] < c.vertex_to_col[_n2])
                                can_n = _n1;
                            else
                                can_n = _n2;
                            const int orig_test_v = translate_back(_test_v);
                            const int orig_n1 = translate_back(can_n);
                            //assert(debug_str_sz == recovery_strings[orig_test_v].size());
                            recovery_strings[orig_n1].push_back(orig_test_v);
                            for (size_t s = 0; s < recovery_strings[orig_test_v].size(); ++s) {
                                recovery_strings[orig_n1].push_back(recovery_strings[orig_test_v][s]);
                            }
                        }
                    }
                }
                i += c.ptn[i] + 1;
            }
 
            worklist_deg0.reset();
            worklist_deg1.reset();
            touched_color_cache.reset();
        }
 
        // in g, walk from 'start' (degree 2) until vertex of not-degree 2 is reached
        // never walks to 'block', if adjacent to 'start'
        // watch out! won't terminate on cycles
        static int walk_cycle(dejavu::sgraph *g, const int start, const int block, dejavu::markset* path_done,
                              dejavu::worklist* path) {
            int current_vertex = start;
            int last_vertex    = block;
 
            int length = 1;
 
            if(path)      path->push_back(block);
            if(path_done) path_done->set(block);
 
            while(true) {
                if (g->d[current_vertex] != 2) {
                    length = -1;
                    break;
                }
                ++length;
 
                if(path)      path->push_back(current_vertex);
                if(path_done) path_done->set(current_vertex);
 
                const int next_vertex1 = g->e[g->v[current_vertex] + 0];
                const int next_vertex2 = g->e[g->v[current_vertex] + 1];
                int next_vertex = next_vertex1;
                if(next_vertex1 == last_vertex) next_vertex = next_vertex2;
                if(next_vertex  == block)       break;
 
                last_vertex    = current_vertex;
                current_vertex = next_vertex;
            }
 
            return length;
        }
 
        // in g, walk from 'start' (degree 2) until vertex of not-degree 2 is reached
        // never walks to 'block', if adjacent to 'start'
        // watch out! won't terminate on cycles
        static std::pair<int, int> walk_to_endpoint(dejavu::sgraph *g, const int start, const int block,
                                                    dejavu::markset* path_done) {
            int current_vertex = start;
            int last_vertex    = block;
 
            while(g->d[current_vertex] == 2) {
                if(path_done)
                    path_done->set(current_vertex);
 
                const int next_vertex1 = g->e[g->v[current_vertex] + 0];
                const int next_vertex2 = g->e[g->v[current_vertex] + 1];
                int next_vertex = next_vertex1;
                if(next_vertex1 == last_vertex)
                    next_vertex = next_vertex2;
 
                last_vertex    = current_vertex;
                current_vertex = next_vertex;
            }
 
            return {current_vertex, last_vertex};
        }
 
        void order_edgelist(dejavu::sgraph *g) {
            int k;
            for(k = 1; k < g->v_size && g->v[k-1] <= g->v[k]; ++k);
            if(k == g->v_size) return; // already ordered
 
            memcpy(edge_scratch.get_array(), g->e, g->e_size*sizeof(int));
            int epos = 0;
            for(int i = 0; i < g->v_size; ++i) {
                const int eptr = g->v[i];
                const int deg  = g->d[i];
                g->v[i] = epos;
                for(int j = eptr; j < eptr + deg; ++j) {
                    g->e[epos] = edge_scratch[j];
                    ++epos;
                }
            }
        }
 
        // in g, walk from 'start' (degree 2) until vertex of not-degree 2 is reached
        // never walks to 'block', if adjacent to 'start'
        // watch out! won't terminate on cycles
        static int walk_to_endpoint_collect_path(dejavu::sgraph *g, const int start, const int block,
                                                 dejavu::worklist* path) {
            int current_vertex = start;
            int last_vertex    = block;
 
            while(g->d[current_vertex] == 2) {
                if(path)
                    path->push_back(current_vertex);
 
                const int next_vertex1 = g->e[g->v[current_vertex] + 0];
                const int next_vertex2 = g->e[g->v[current_vertex] + 1];
                int next_vertex = next_vertex1;
                if(next_vertex1 == last_vertex)
                    next_vertex = next_vertex2;
 
                last_vertex    = current_vertex;
                current_vertex = next_vertex;
            }
 
            dej_assert(g->d[current_vertex] != 2);
            return current_vertex;
        }
 
        // color-wise degree2 unique endpoint algorithm
        void red_deg2_unique_endpoint(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1 || h_deact_deg2)
                return;
 
            dejavu::markset color_test(g->v_size);
            dejavu::markset color_unique(g->v_size);
 
            //coloring col;
            g->initialize_coloring(&c, colmap);
 
            add_edge_buff_act.reset();
            for (int i = 0; i < g->v_size; ++i) {
                add_edge_buff[i].clear();
            }
 
            del.reset();
 
            worklist_deg1.reset();
 
            dejavu::worklist endpoint_cnt(g->v_size);
            for (int i = 0; i < g->v_size; ++i) {
                endpoint_cnt.push_back(0);
            }
 
            dejavu::markset  path_done(g->v_size);
            dejavu::worklist color_pos(g->v_size);
            dejavu::worklist filter(g->v_size);
            dejavu::worklist path_list(g->v_size);
            dejavu::worklist path(g->v_size);
            dejavu::worklist connected_paths(g->e_size);
            dejavu::worklist connected_endpoints(g->e_size);
 
            for (int i = 0; i < g->v_size; ++i) {
                if(g->d[i] == 2) {
                    if(path_done.get(i))
                        continue;
                    const int n1 = g->e[g->v[i] + 0];
                    const int n2 = g->e[g->v[i] + 1];
                    // n1 or n2 is endpoint? then walk to other endpoint and collect information...
                    if(g->d[n1] != 2) {
                        const auto other_endpoint = walk_to_endpoint(g, i, n1, &path_done);
                        if(other_endpoint.first == n1) // big self-loop
                            continue;
                        connected_paths[g->v[n1] + endpoint_cnt[n1]]     = i;
                        connected_endpoints[g->v[n1] + endpoint_cnt[n1]] = other_endpoint.first;
                        ++endpoint_cnt[n1];
                        connected_paths[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] =
                                other_endpoint.second;
                        connected_endpoints[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] = n1;
                        ++endpoint_cnt[other_endpoint.first];
                        dej_assert(other_endpoint.first != n1);
                        dej_assert(g->d[other_endpoint.first] != 2);
                        dej_assert(n1 != other_endpoint.first);
                    } else if(g->d[n2] != 2) {
                        const auto other_endpoint = walk_to_endpoint(g, i, n2, &path_done);
                        if(other_endpoint.first == n2) // big self-loop
                            continue;
                        connected_paths[g->v[n2] + endpoint_cnt[n2]] = i;
                        connected_endpoints[g->v[n2] + endpoint_cnt[n2]] = other_endpoint.first;
                        ++endpoint_cnt[n2];
                        connected_paths[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] =
                                other_endpoint.second;
                        connected_endpoints[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] = n2;
                        ++endpoint_cnt[other_endpoint.first];
                        dej_assert(other_endpoint.first != n2);
                        dej_assert(g->d[other_endpoint.first] != 2);
                        dej_assert(n2 != other_endpoint.first);
                    }
                }
            }
 
            path_done.reset();
 
            // iterate over color classes
            for(int i = 0; i < g->v_size;) {
                const int color      = i;
                const int color_size = c.ptn[color] + 1;
                i += color_size;
                const int test_vertex = c.lab[color];
                const int endpoints   = endpoint_cnt[test_vertex];
 
                for(int j = 0; j < color_size; ++j) {
                    dej_assert(endpoint_cnt[c.lab[color + j]] == endpoints);
                }
 
                if(endpoints == 0) continue;
 
 
                color_test.reset();
                color_unique.reset();
 
                // check which neighbouring paths have unique color
                for(int j = 0; j < endpoints; ++j) {
                    const int neighbour     = connected_paths[g->v[test_vertex] + j];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if(color_test.get(neighbour_col)) {
                        color_unique.set(neighbour_col); // means "not unique"
                    }
                    color_test.set(neighbour_col);
                }
 
                filter.reset();
                // filter to indices with unique colors
                for(int j = 0; j < endpoints; ++j) {
                    const int neighbour     = connected_paths[g->v[test_vertex] + j];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if(!color_unique.get(neighbour_col)) { // if unique
                        filter.push_back(j); // add index to filter
                    }
                }
 
                path.reset();
                color_test.reset();
                color_unique.reset();
 
                // check which neighbouring endpoints have unique color and are NOT connected to `color`
                for(int j = 0; j < g->d[test_vertex]; ++j) {
                    const int neighbour     = g->e[g->v[test_vertex] + j];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    color_test.set(neighbour_col);
                }
 
                // also self-connections of color are forbidden
                color_test.set(color);
 
                for(int k = 0; k < filter.cur_pos; ++k) {
                    const int j = filter[k];
                    const int neighbour_endpoint = connected_endpoints[g->v[test_vertex] + j];
                    const int neighbour_col      = c.vertex_to_col[neighbour_endpoint];
                    if(color_test.get(neighbour_col)) {
                        color_unique.set(neighbour_col);
                    }
                    color_test.set(neighbour_col);
                }
 
                // filter to indices with unique colors
                int write_pos = 0;
                for(int k = 0; k < filter.cur_pos; ++k) {
                    const int j = filter[k];
                    const int neighbour     = connected_endpoints[g->v[test_vertex] + j];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if(!color_unique.get(neighbour_col)) {
                        filter[write_pos] = j;
                        ++write_pos;
                    }
                }
                filter.cur_pos = write_pos;
 
 
                // encode filter into color_unique, such that we can detect corresponding neighbours for all vertices
                // of color class
                color_unique.reset();
                for(int k = 0; k < filter.cur_pos; ++k) {
                    const int j = filter[k];
                    const int filter_to_col = c.vertex_to_col[connected_paths[g->v[test_vertex] + j]];
                    color_unique.set(filter_to_col);
                    filter[k] = filter_to_col;
                    //color_unique.set(col.vertex_to_col[connected_endpoints[filter[k]]]);
                }
 
                // order colors here already to access color_pos in O(1) later
                filter.sort();
                for(int k = 0; k < filter.cur_pos; ++k) {
                    color_pos[filter[k]] = k;
                }
 
                const int num_paths = filter.cur_pos;
 
                #if defined(DEJDEBUG) &&  !defined(NDEBUG)
                int reduced_verts = 0;
                int reduced_verts_last = 0;
                #endif
 
                // do next steps for all vertices in color classes...
                for(int j = 0; j < color_size; ++j) {
                    const int vertex = c.lab[color + j];
 
                    int sanity_check = 0;
                    // paths left in filter (color_unique) are now collected for reduction
                    color_test.reset();
                    for(int k = 0; k < g->d[vertex]; ++k) {
                        const int neighbour     = g->e[g->v[vertex] + k];
                        const int neighbour_col = c.vertex_to_col[neighbour];
                        if(color_unique.get(neighbour_col)) {
                            const int pos = color_pos[neighbour_col];
                            path_list[pos] = neighbour;
                            ++sanity_check;
                            dej_assert(!color_test.get(neighbour_col));
                            color_test.set(neighbour_col);
                        }
                    }
 
                    dej_assert(sanity_check == num_paths);
 
                    for(int k = 0; k < num_paths; ++k) {
                        const int path_start_vertex = path_list[k];
                        dej_assert(g->d[path_start_vertex] == 2);
                        if(path_done.get(path_start_vertex)) {
                            continue;
                        }
 
                        // get path and endpoint
                        path.reset();
                        const int other_endpoint =
                                walk_to_endpoint_collect_path(g, path_start_vertex, vertex, &path);
                        dej_assert(path.cur_pos > 0);
                        dej_assert(vertex != other_endpoint);
                        dej_assert(other_endpoint != path_start_vertex);
 
                        // mark path for deletion
                        for(int l = 0; l < path.cur_pos; ++l) {
                            const int del_v = path[l];
                            dej_assert(g->d[del_v] == 2);
                            dej_assert(!del.get(del_v));
                            del.set(del_v);
                            dej_assert(!path_done.get(del_v));
                            path_done.set(del_v);
                        }
 
                        // connect endpoints of path with new edge
                        dej_assert(!is_adjacent(g, vertex, other_endpoint));
                        add_edge_buff[other_endpoint].push_back(vertex);
                        add_edge_buff_act.set(other_endpoint);
                        add_edge_buff[vertex].push_back(other_endpoint);
                        add_edge_buff_act.set(vertex);
 
                        // write path into recovery_strings
                        const int unique_endpoint_orig = translate_back(vertex);
                        // attach all represented vertices of path to unique_endpoint_orig in canonical fashion
                        // path.sort_after_map(colmap); // should not be necessary
                        recovery_strings[unique_endpoint_orig].reserve(
                                2*(recovery_strings[unique_endpoint_orig].size() + path.cur_pos));
                        for (int l = 0; l < path.cur_pos; ++l) {
                            dej_assert(path[l] >= 0);
                            dej_assert(path[l] < g->v_size);
                            const int path_v_orig = translate_back(path[l]);
                            dej_assert(path_v_orig >= 0);
                            dej_assert(path_v_orig < domain_size);
                            recovery_strings[unique_endpoint_orig].push_back(path_v_orig);
                            recovery_strings[unique_endpoint_orig].insert(
                                    recovery_strings[unique_endpoint_orig].end(),
                                    recovery_strings[path_v_orig].begin(),
                                    recovery_strings[path_v_orig].end());
                        }
                    }
 
                    #if defined(DEJDEBUG) &&  !defined(NDEBUG)
                    reduced_verts = static_cast<int>(recovery_strings[translate_back(vertex)].size());
                    if(j > 0) {
                        dej_assert(reduced_verts == reduced_verts_last);
                    }
                    reduced_verts_last = reduced_verts;
                    #endif
                }
            }
        }
 
        void recovery_attach_to_edge(int v1, int v2, std::vector<int>& recovery) {
            const int pos = static_cast<int>(recovery_edge_attached.size());
            dej_assert(recovery.size() != 0);
            const int len = static_cast<int>(recovery.size());
            ++edge_attached_information;
            if(len == 1) {
                // special code for stored path of length 1 (does not use recovery_edge_attached)
                recovery_edge_adjacent[v1].emplace_back(v2, recovery[0], len);
                recovery_edge_adjacent[v2].emplace_back(v1, recovery[0], len);
            } else {
                // general-purpose code for longer paths
                for (int i = 0; i < len; ++i) { recovery_edge_attached.push_back(recovery[i]); }
                recovery_edge_adjacent[v1].emplace_back(v2, pos, len);
                recovery_edge_adjacent[v2].emplace_back(v1, pos, len);
            }
        }
 
        void recovery_attached_edges_to_automorphism(int v, int* aut, dejavu::ds::worklist* aut_supp,
                                                     dejavu::ds::worklist* help_array1,
                                                     dejavu::ds::worklist* help_array2) {
            if(recovery_edge_adjacent[v].empty()) return;
            const int v_map    = aut[v];
 
            for(auto & j : recovery_edge_adjacent[v]) {
                const int other = std::get<0>(j);
                const int id    = std::get<1>(j);
 
                dej_assert(other != v);
                const int other_map = aut[other];
 
                (*help_array2)[other_map] = id;
            }
 
            dej_assert(recovery_edge_adjacent[v].size() == recovery_edge_adjacent[v_map].size());
 
            for(auto & j : recovery_edge_adjacent[v_map]) {
                const int other_map = std::get<0>(j);
                const int id        = std::get<1>(j);
                const int len       = std::get<2>(j);
 
                dej_assert(other_map != v_map);
                const int orig_id = (*help_array2)[other_map];
                if(id == orig_id) continue;
 
                // endpoints of the mapped edges:
                // e1: v     -> v_map
                // e2: other -> other_map
                // aut maps e1 -> e2
 
                if(len == 1) {
                    // special code for stored path of length 1 (does not use recovery_edge_attached)
                    const int v_from = orig_id;
                    const int v_to   = id;
                    if (v_from != v_to && aut[v_from] == v_from) {
                        aut[v_from] = v_to;
                        aut_supp->push_back(v_from);
                    }
                } else {
                    // general-purpose code for longer paths
                    for (int k = 0; k < len; ++k) {
                        const int v_from = recovery_edge_attached[orig_id + k];
                        const int v_to   = recovery_edge_attached[id + k];
                        if (v_from != v_to && aut[v_from] == v_from) {
                            aut[v_from] = v_to;
                            aut_supp->push_back(v_from);
                        }
                    }
                }
            }
        }
 
        void red_deg2_densify(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1 || h_deact_deg2) return;
 
            dejavu::markset color_test(g->v_size);
            dejavu::markset color_unique(g->v_size);
 
            //coloring col;
            g->initialize_coloring(&c, colmap);
            add_edge_buff_act.reset();
            for (int i = 0; i < g->v_size; ++i) {
                add_edge_buff[i].clear();
            }
 
            del.reset();
            worklist_deg1.reset();
 
            dejavu::worklist endpoint_cnt(g->v_size);
            for (int i = 0; i < g->v_size; ++i) endpoint_cnt.push_back(0);
 
            dejavu::markset  path_done(g->v_size);
            dejavu::worklist color_pos(g->v_size);
            dejavu::worklist color_deg(g->v_size);
            dejavu::worklist hub_vertex_position(g->v_size);
            dejavu::worklist color_canonical_v(g->v_size);
            dejavu::worklist filter(g->v_size);
            dejavu::worklist path_list(g->v_size);
            dejavu::worklist path(g->v_size);
            dejavu::worklist connected_paths(g->e_size);
            dejavu::worklist connected_endpoints(g->e_size);
            dejavu::markset  duplicate_endpoints(g->v_size);
 
            // collect and count endpoints
            int total_paths = 0;
 
            for (int i = 0; i < g->v_size; ++i) {
                if(g->d[i] == 2) {
                    hub_vertex_position[i] = -1;
                    if(!recovery_strings[translate_back(i)].empty()) {
                        path_done.set(i); // not ideal...
                    }
                    if(path_done.get(i))
                        continue;
                    const int n1 = g->e[g->v[i] + 0];
                    const int n2 = g->e[g->v[i] + 1];
                    // n1 or n2 is endpoint? then walk to other endpoint and collect information...
                    if(g->d[n1] != 2) {
                        const auto other_endpoint = walk_to_endpoint(g, i, n1, &path_done);
                        if(other_endpoint.first == n1) // big self-loop
                            continue;
                        connected_paths[g->v[n1] + endpoint_cnt[n1]]     = i;
                        connected_endpoints[g->v[n1] + endpoint_cnt[n1]] = other_endpoint.first;
                        ++endpoint_cnt[n1];
                        connected_paths[g->v[other_endpoint.first]     + endpoint_cnt[other_endpoint.first]] =
                                other_endpoint.second;
                        connected_endpoints[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] = n1;
                        ++endpoint_cnt[other_endpoint.first];
                        dej_assert(other_endpoint.first != n1);
                        dej_assert(g->d[other_endpoint.first] != 2);
                        dej_assert(n1 != other_endpoint.first);
                        ++total_paths;
                    } else if(g->d[n2] != 2) {
                        const auto other_endpoint = walk_to_endpoint(g, i, n2, &path_done);
                        if(other_endpoint.first == n2) // big self-loop
                            continue;
                        connected_paths[g->v[n2] + endpoint_cnt[n2]] = i;
                        connected_endpoints[g->v[n2] + endpoint_cnt[n2]] = other_endpoint.first;
                        ++endpoint_cnt[n2];
                        connected_paths[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] =
                                other_endpoint.second;
                        connected_endpoints[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] = n2;
                        ++endpoint_cnt[other_endpoint.first];
                        dej_assert(other_endpoint.first != n2);
                        dej_assert(g->d[other_endpoint.first] != 2);
                        dej_assert(n2 != other_endpoint.first);
                        ++total_paths;
                    }
                }
            }
 
            // let's not apply this reduction unless there is a substantial amount of degree 2 vertices -- performing
            // the technique adds some bloat to the lifting process, so we should not overdo it
            if(total_paths < g->v_size/10) return;
 
            recovery_edge_adjacent.resize(domain_size);
            path_done.reset();
 
            // iterate over color classes
            for(int i = 0; i < g->v_size;) {
                const int color      = i;
                const int color_size = c.ptn[color] + 1;
                i += color_size;
                const int test_vertex = c.lab[color];
                const int endpoints   = endpoint_cnt[test_vertex];
 
                for(int j = 0; j < color_size; ++j) {
                    dej_assert(endpoint_cnt[c.lab[color + j]] == endpoints);
                }
 
                // only want to look at endpoints of paths
                if(endpoints == 0) continue;
 
                color_test.reset();
                color_unique.reset();
 
                // only interested in paths with non-unique color, otherwise we could use unique endpoint algorithm
                // check which neighbouring paths have non-unique color
                for(int j = 0; j < endpoints; ++j) {
                    const int neighbour     = connected_paths[g->v[test_vertex] + j];
 
                    dej_assert(g->d[neighbour] == 2);
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if(color_test.get(neighbour_col) && !path_done.get(neighbour)) {
                        ++color_deg[neighbour_col];
                        color_unique.set(neighbour_col); // means "not unique"
                    } else {
                        color_deg[neighbour_col] = 1;
                    }
                    color_test.set(neighbour_col);
                }
 
                // Make sure there are no duplicate endpoints, i.e., paths going from the same vertex to the same
                // vertex!
                bool duplicate_endpoint_flag = false;
                for(int j = 0; j < color_size; ++j) {
                    const int endpt_test_vertex = c.lab[color + j];
                    duplicate_endpoints.reset();
                    duplicate_endpoints.set(endpt_test_vertex); // no self-loops! tricky case
                    dej_assert(endpoint_cnt[endpt_test_vertex] == endpoints);
                    for (int k = 0; k < endpoints; ++k) {
                        const int endpoint = connected_endpoints[g->v[endpt_test_vertex] + k];
                        if(duplicate_endpoints.get(endpoint)) {
                            duplicate_endpoint_flag = true;
                            break;
                        }
                        duplicate_endpoints.set(endpoint);
                    }
                    if(duplicate_endpoint_flag) break;
                }
 
                // There were duplicate endpoints, so we should stop with this color class!
                if(duplicate_endpoint_flag) continue;
 
                filter.reset();
                // Filter to indices with non-unique colors
                for(int j = 0; j < endpoints; ++j) {
                    const int neighbour     = connected_paths[g->v[test_vertex] + j];
                    dej_assert(g->d[neighbour] == 2);
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    const int endpoint     = connected_endpoints[g->v[test_vertex] + j];
                    dej_assert(g->d[endpoint] != 2);
                    const int endpoint_col = c.vertex_to_col[endpoint];
                    const int endpoint_col_sz = c.ptn[endpoint_col] + 1;
                    if(color_unique.get(neighbour_col) && color != endpoint_col &&
                       endpoint_col_sz >= color_size) { // if unique
                        filter.push_back(j); // add index to filter
                    }
                }
 
 
                path.reset();
                color_test.reset();
                color_unique.reset();
 
                int use_pos = 0;
                for(int k = 0; k < filter.cur_pos; ++k) {
                    const int filter_col = c.vertex_to_col[connected_paths[g->v[test_vertex] + filter[k]]];
                    if(!color_unique.get(filter_col)) {
                        color_pos[filter_col] = use_pos;
                        color_unique.set(filter_col);
                        use_pos += color_deg[filter_col];
                    }
                }
 
                const int num_paths = filter.cur_pos;
 
                // do next steps for all vertices in color classes...
                for(int j = 0; j < color_size; ++j) {
                    const int vertex = c.lab[color + j];
 
                    // reset color_deg to 0
                    for(int k = 0; k < filter.cur_pos; ++k) {
                        const int filter_col = c.vertex_to_col[connected_paths[g->v[test_vertex] + filter[k]]];
                        color_deg[filter_col] = 0;
                    }
 
                    int sanity_check_num_paths = 0;
 
                    for(int k = 0; k < g->d[vertex]; ++k) {
                        const int neighbour     = g->e[g->v[vertex] + k];
                        const int neighbour_col = c.vertex_to_col[neighbour];
                        if(color_unique.get(neighbour_col)) {
                            const int pos = color_pos[neighbour_col] + color_deg[neighbour_col];
                            path_list[pos] = neighbour;
                            dej_assert(g->d[neighbour] == 2);
                            color_canonical_v[neighbour_col] = neighbour; // "hub vertex" used to emulate all the paths
                            ++color_deg[neighbour_col];
                            ++sanity_check_num_paths;
                        }
                    }
 
                    dej_assert(num_paths == sanity_check_num_paths);
 
                    for(int k = 0; k < num_paths; ++k) {
                        const int path_start_vertex = path_list[k];
                        dej_assert(g->d[path_start_vertex] == 2);
                        if(path_done.get(path_start_vertex)) {
                            continue;
                        }
 
                        // get path and endpoint
                        path.reset();
                        const int other_endpoint = walk_to_endpoint_collect_path(g, path_start_vertex, vertex, &path);
                        dej_assert(path.cur_pos > 0);
                        dej_assert(vertex != other_endpoint);
                        dej_assert(other_endpoint != path_start_vertex);
 
 
                        const int hub_vertex = color_canonical_v[c.vertex_to_col[path_start_vertex]];
                        dej_assert(hub_vertex >= 0);
                        dej_assert(!del.get(hub_vertex));
                        dej_assert(g->d[hub_vertex] == 2);
                        dej_assert(is_adjacent(g, hub_vertex, vertex));
 
                        // mark path for deletion
                        for(int l = 0; l < path.cur_pos; ++l) {
                            const int del_v = path[l];
                            dej_assert(g->d[del_v] == 2);
                            dej_assert(!del.get(del_v));
                            del.set(del_v);
                            dej_assert(!path_done.get(del_v));
                            path_done.set(del_v);
                        }
                        del.unset(hub_vertex); // don't delete canonical vertex
                        dej_assert(!del.get(hub_vertex));
 
                        // connect endpoints of path to hub vertex (the canonical vertex)
                        // make sure not do double up on hub_vertex original path
                        if(!is_adjacent(g, hub_vertex, other_endpoint)) {
                            add_edge_buff[other_endpoint].push_back(hub_vertex);
                            add_edge_buff_act.set(other_endpoint);
                            add_edge_buff[hub_vertex].push_back(other_endpoint);
                            add_edge_buff_act.set(hub_vertex);
                        }
 
                        // write path into recovery_strings
                        const int hub_vertex_orig = translate_back(hub_vertex);
 
                        // mark hub vertex in recovery string such that it gets mapped to nothing
                        // TODO should this be handled by translate_back? we can "translate_back" to \epsilon?
                        // TODO just collect hub vertices and "handle this at the end"?
                        if(recovery_strings[hub_vertex_orig].empty())
                            recovery_strings[hub_vertex_orig].push_back(INT32_MAX);
 
                        std::vector<int> recovery;
                        for (int l = 0; l < path.cur_pos; ++l) {
                            dej_assert(path[l] >= 0);
                            dej_assert(path[l] < g->v_size);
                            const int path_v_orig = translate_back(path[l]);
                            //if(hub_vertex_orig == path_v_orig) continue;
                            dej_assert(path_v_orig >= 0);
                            dej_assert(path_v_orig < domain_size);
                            recovery.push_back(path_v_orig);
                            dej_assert(path_v_orig != hub_vertex_orig? recovery_strings[path_v_orig].size() == 0: true);
                        }
 
                        recovery_attach_to_edge(translate_back(vertex), translate_back(other_endpoint), recovery);
                    }
                }
            }
        }
 
        void red_deg2_densifysc1(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1 || h_deact_deg2) return;
 
            //coloring col;
            g->initialize_coloring(&c, colmap);
            add_edge_buff_act.reset();
            for (int i = 0; i < g->v_size; ++i) add_edge_buff[i].clear();
 
            del.reset();
            worklist_deg1.reset();
 
            dejavu::worklist endpoint_cnt(g->v_size);
            for (int i = 0; i < g->v_size; ++i) endpoint_cnt.push_back(0);
 
            dejavu::markset  path_done(g->v_size);
            dejavu::worklist filter(g->v_size);
            dejavu::worklist path_list(g->v_size);
            dejavu::worklist path(g->v_size);
            dejavu::worklist connected_paths(g->e_size);
            dejavu::worklist connected_endpoints(g->e_size);
            dejavu::markset  duplicate_endpoints(g->v_size);
 
            // collect and count endpoints
            int total_paths = 0;
 
            for (int i = 0; i < g->v_size; ++i) {
                if(g->d[i] == 2) {
                    if(path_done.get(i))
                        continue;
                    const int n1 = g->e[g->v[i] + 0];
                    const int n2 = g->e[g->v[i] + 1];
                    // n1 or n2 is endpoint? then walk to other endpoint and collect information...
                    if(g->d[n1] != 2 && g->d[n2] != 2) { // only want paths of length 1
                        const auto other_endpoint = walk_to_endpoint(g, i, n1, &path_done);
                        if(other_endpoint.first == n1) // big self-loop
                            continue;
                        connected_paths[g->v[n1] + endpoint_cnt[n1]]     = i;
                        connected_endpoints[g->v[n1] + endpoint_cnt[n1]] = other_endpoint.first;
                        ++endpoint_cnt[n1];
                        connected_paths[g->v[other_endpoint.first]     + endpoint_cnt[other_endpoint.first]] =
                                other_endpoint.second;
                        connected_endpoints[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] = n1;
                        ++endpoint_cnt[other_endpoint.first];
                        dej_assert(other_endpoint.first != n1);
                        dej_assert(g->d[other_endpoint.first] != 2);
                        dej_assert(n1 != other_endpoint.first);
                        ++total_paths;
                    }
                }
            }
 
            // let's not apply this reduction unless there is a substantial amount of degree 2 vertices -- performing
            // the technique adds some bloat to the lifting process, so we should not overdo it
            if(total_paths < g->v_size/10) return;
 
            recovery_edge_adjacent.resize(domain_size);
            path_done.reset();
 
            // iterate over color classes
            for(int i = 0; i < g->v_size;) {
                const int color      = i;
                const int color_size = c.ptn[color] + 1;
                i += color_size;
                const int test_vertex = c.lab[color];
                const int endpoints   = endpoint_cnt[test_vertex];
 
                for(int j = 0; j < color_size; ++j) {
                    dej_assert(endpoint_cnt[c.lab[color + j]] == endpoints);
                }
 
                // only want to look at endpoints of paths
                if(endpoints == 0) continue;
 
                // only interested in paths with non-unique color, otherwise we could use unique endpoint algorithm
                // check which neighbouring paths have non-unique color
                int color_deg_single = 0;
                bool only_one = true;
                int relevant_neighbour_color = -1;
                for(int j = 0; j < endpoints; ++j) {
                    const int neighbour     = connected_paths[g->v[test_vertex] + j];
                    const int endpoint      = connected_endpoints[g->v[test_vertex] + j];
                    dej_assert(g->d[endpoint] != 2);
                    dej_assert(g->d[neighbour] == 2);
                    const int endpoint_col = c.vertex_to_col[endpoint];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if(!path_done.get(neighbour) && endpoint_col == color) {
                        if(relevant_neighbour_color == -1) {
                            relevant_neighbour_color = neighbour_col;
                        }
                        only_one = only_one && (relevant_neighbour_color == neighbour_col);
                        color_deg_single += 1;
                    }
                }
 
                if(!only_one || color_deg_single == 0) continue;
 
                bool duplicate_endpoint_flag = false;
                // Make sure there is no self-connection to own color already!
                for(int j = g->v[test_vertex]; j < g->v[test_vertex] + g->d[test_vertex]; ++j) {
                    const int neighbour = g->e[j];
                    if(c.vertex_to_col[neighbour] == color) {
                        duplicate_endpoint_flag = true;
                        break;
                    }
                }
 
                // Make sure there are no duplicate endpoints, i.e., paths going from the same vertex to the same
                // vertex!
                for(int j = 0; j < color_size; ++j) {
                    const int endpt_test_vertex = c.lab[color + j];
                    duplicate_endpoints.reset();
                    duplicate_endpoints.set(endpt_test_vertex); // no self-loops! tricky case
                    dej_assert(endpoint_cnt[endpt_test_vertex] == endpoints);
                    for (int k = 0; k < endpoints; ++k) {
                        const int neighbour = connected_paths[g->v[endpt_test_vertex] + k];
                        const int neighbour_col = c.vertex_to_col[neighbour];
                        if(neighbour_col != relevant_neighbour_color) continue;
                        dej_assert(g->d[neighbour] == 2);
                        const int endpoint = connected_endpoints[g->v[endpt_test_vertex] + k];
                        dej_assert(c.vertex_to_col[endpoint] == color);
                        if(duplicate_endpoints.get(endpoint)) {
                            duplicate_endpoint_flag = true;
                            break;
                        }
                        duplicate_endpoints.set(endpoint);
                    }
                    if(duplicate_endpoint_flag) break;
                }
 
                // There were duplicate endpoints, so we should stop with this color class!
                if(duplicate_endpoint_flag) continue;
 
                filter.reset();
                // Filter to indices with non-unique colors
                for(int j = 0; j < endpoints; ++j) {
                    const int neighbour     = connected_paths[g->v[test_vertex] + j];
                    dej_assert(g->d[neighbour] == 2);
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if(neighbour_col == relevant_neighbour_color) { // if unique
                        dej_assert(neighbour_col == relevant_neighbour_color);
                        filter.push_back(j); // add index to filter
                    }
                }
 
                path.reset();
                const int num_paths = filter.cur_pos;
 
                // do next steps for all vertices in color classes...
                for(int j = 0; j < color_size; ++j) {
                    const int vertex = c.lab[color + j];
 
                    // reset color_deg to 0
                    int next_pos = 0;
                    int sanity_check_num_paths = 0;
 
                    for(int k = 0; k < g->d[vertex]; ++k) {
                        const int neighbour     = g->e[g->v[vertex] + k];
                        const int neighbour_col = c.vertex_to_col[neighbour];
                        if(neighbour_col == relevant_neighbour_color) {
                            const int pos = next_pos;
                            path_list[pos] = neighbour;
                            dej_assert(g->d[neighbour] == 2);
                            ++next_pos;
                            ++sanity_check_num_paths;
                        }
                    }
 
                    dej_assert(num_paths == sanity_check_num_paths);
 
                    for(int k = 0; k < num_paths; ++k) {
                        const int path_start_vertex = path_list[k];
                        dej_assert(g->d[path_start_vertex] == 2);
                        if(path_done.get(path_start_vertex)) {
                            continue;
                        }
 
                        // get path and endpoint
                        path.reset();
                        const int other_endpoint = walk_to_endpoint_collect_path(g, path_start_vertex, vertex, &path);
                        dej_assert(path.cur_pos > 0);
                        dej_assert(vertex != other_endpoint);
                        dej_assert(other_endpoint != path_start_vertex);
 
                        // mark path for deletion
                        for(int l = 0; l < path.cur_pos; ++l) {
                            const int del_v = path[l];
                            dej_assert(g->d[del_v] == 2);
                            dej_assert(!del.get(del_v));
                            del.set(del_v);
                            dej_assert(!path_done.get(del_v));
                            path_done.set(del_v);
                        }
 
                        // connect endpoints of path to hub vertex (the canonical vertex)
                        // make sure not do double up on hub_vertex original path
                        if(!is_adjacent(g, vertex, other_endpoint)) {
                            add_edge_buff[other_endpoint].push_back(vertex);
                            add_edge_buff_act.set(other_endpoint);
                            add_edge_buff[vertex].push_back(other_endpoint);
                            add_edge_buff_act.set(vertex);
                        }
 
                        // write path into recovery_strings
                        std::vector<int> recovery;
                        for (int l = 0; l < path.cur_pos; ++l) {
                            dej_assert(path[l] >= 0);
                            dej_assert(path[l] < g->v_size);
                            const int path_v_orig = translate_back(path[l]);
                            dej_assert(path_v_orig >= 0);
                            dej_assert(path_v_orig < domain_size);
                            recovery.push_back(path_v_orig);
                            for(int f = 0; f < static_cast<int>(recovery_strings[path_v_orig].size()); ++f) {
                                recovery.push_back(recovery_strings[path_v_orig][f]);
                            }
                        }
 
                        recovery_attach_to_edge(translate_back(vertex), translate_back(other_endpoint), recovery);
                    }
                }
            }
        }
 
        // color cycles according to their size
        // remove uniform colored cycles
        bool red_deg2_color_cycles(dejavu::sgraph *g, int *colmap) {
            //coloring col;
            g->initialize_coloring(&c, colmap);
            for(int i = 0; i < g->v_size; ++i) {
                colmap[i] = c.vertex_to_col[i];
            }
 
            int cycles_recolored = 0;
 
            dejavu::markset  path_done(g->v_size);
            dejavu::worklist path(g->v_size);
 
            for (int i = 0; i < g->v_size; ++i) {
                if(g->d[i] == 2) {
                    if(path_done.get(i))
                        continue;
                    const int n1 = g->e[g->v[i] + 0];
                    const int n2 = g->e[g->v[i] + 1];
 
                    if(g->d[n1] != 2 || g->d[n2] != 2)
                        continue;
 
                    path.reset();
                    const int cycle_length = walk_cycle(g, i, n1, &path_done, &path);
                    if(cycle_length == -1) {
                        continue;
                    } else {
                        cycles_recolored += 1;
                        for(int j = 0; j < path.cur_pos; ++j) {
                            colmap[path[j]] = colmap[path[j]] + g->v_size * cycle_length; // maybe this should receive a more elegant solution
                        }
                    }
                }
            }
 
            if(cycles_recolored > 0) {
                g->initialize_coloring(&c, colmap);
                bool has_discrete = false;
                for(int i = 0; i < g->v_size && !has_discrete; ) {
                    const int col_sz = c.ptn[i] + 1;
                    has_discrete = has_discrete || (col_sz == 1);
                    i += col_sz;
                }
                for(int i = 0; i < g->v_size; ++i) {
                    colmap[i] = c.vertex_to_col[i];
                }
                return has_discrete;
            }
            return false;
        }
 
        void red_deg2_trivial_connect(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1 || h_deact_deg2)
                return;
 
            dejavu::ds::markset color_test(g->v_size);
 
            dejavu::ds::markset color_unique(g->v_size);
 
            //coloring col;
            g->initialize_coloring(&c, colmap);
 
            add_edge_buff_act.reset();
            for (int i = 0; i < g->v_size; ++i) add_edge_buff[i].clear();
 
            del.reset();
 
            worklist_deg1.reset();
 
            dejavu::ds::worklist endpoint_cnt(g->v_size);
            for (int i = 0; i < g->v_size; ++i) endpoint_cnt.push_back(0);
 
            dejavu::ds::markset path_done(g->v_size);
            dejavu::ds::worklist color_pos(g->v_size);
            dejavu::ds::worklist not_unique(2*g->v_size);
            dejavu::ds::worklist not_unique_analysis(g->v_size);
            dejavu::ds::worklist path_list(g->v_size);
            dejavu::ds::worklist path(g->v_size);
            dejavu::ds::worklist connected_paths(g->e_size);
            dejavu::ds::worklist connected_endpoints(g->e_size);
            dejavu::ds::worklist neighbour_list(g->v_size);
            dejavu::ds::worklist neighbour_to_endpoint(g->v_size);
 
            for (int i = 0; i < g->v_size; ++i) {
                if(g->d[i] == 2) {
                    if(path_done.get(i))
                        continue;
                    const int n1 = g->e[g->v[i] + 0];
                    const int n2 = g->e[g->v[i] + 1];
                    // n1 or n2 is endpoint? then walk to other endpoint and collect information...
                    if(g->d[n1] != 2) {
                        const auto other_endpoint = walk_to_endpoint(g, i, n1, &path_done);
                        if(other_endpoint.first == n1) // big self-loop
                            continue;
                        connected_paths[g->v[n1] + endpoint_cnt[n1]]     = i;
                        connected_endpoints[g->v[n1] + endpoint_cnt[n1]] = other_endpoint.first;
                        ++endpoint_cnt[n1];
                        connected_paths[g->v[other_endpoint.first]     + endpoint_cnt[other_endpoint.first]] =
                                other_endpoint.second;
                        connected_endpoints[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] = n1;
                        ++endpoint_cnt[other_endpoint.first];
                        dej_assert(other_endpoint.first != n1);
                        dej_assert(g->d[other_endpoint.first] != 2);
                        dej_assert(n1 != other_endpoint.first);
                    } else if(g->d[n2] != 2) {
                        const auto other_endpoint = walk_to_endpoint(g, i, n2, &path_done);
                        if(other_endpoint.first == n2) // big self-loop
                            continue;
                        connected_paths[g->v[n2] + endpoint_cnt[n2]] = i;
                        connected_endpoints[g->v[n2] + endpoint_cnt[n2]] = other_endpoint.first;
                        ++endpoint_cnt[n2];
                        connected_paths[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] =
                                other_endpoint.second;
                        connected_endpoints[g->v[other_endpoint.first] + endpoint_cnt[other_endpoint.first]] = n2;
                        ++endpoint_cnt[other_endpoint.first];
                        dej_assert(other_endpoint.first != n2);
                        dej_assert(g->d[other_endpoint.first] != 2);
                        dej_assert(n2 != other_endpoint.first);
                    }
                }
            }
 
            path_done.reset();
 
            // iterate over color classes
            for(int i = 0; i < g->v_size;) {
                const int color = i;
                const int color_size = c.ptn[color] + 1;
                i += color_size;
                const int test_vertex = c.lab[color];
                const int endpoints = endpoint_cnt[test_vertex];
 
                if (endpoints == 0) {
                    continue;
                }
 
                color_test.reset();
                color_unique.reset();
 
                // check which neighbouring paths have unique color
                for (int j = 0; j < endpoints; ++j) {
                    const int neighbour = connected_paths[g->v[test_vertex] + j];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if (color_test.get(neighbour_col)) {
                        color_unique.set(neighbour_col); // means "not unique"
                    }
                    color_test.set(neighbour_col);
                }
 
                not_unique.reset();
                // filter to indices with unique colors
                for (int j = 0; j < endpoints; ++j) {
                    const int neighbour = connected_paths[g->v[test_vertex] + j];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    if (color_unique.get(neighbour_col)) { // if not unique
                        not_unique.push_back(neighbour);
                        dej_assert(connected_endpoints[g->v[test_vertex] + j] >= 0);
                        dej_assert(connected_endpoints[g->v[test_vertex] + j] < g->v_size);
                        not_unique.push_back(connected_endpoints[g->v[test_vertex] + j]);
                    }
                }
 
                color_test.reset();
                color_unique.reset();
 
                // remove trivial connections
                for (int kk = 0; kk < not_unique.cur_pos; kk += 2) {
                    const int endpoint = not_unique[kk + 1];
                    const int endpoint_col = c.vertex_to_col[endpoint];
                    const int neighbour = not_unique[kk];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    not_unique_analysis[endpoint_col]  = 0;
                    not_unique_analysis[neighbour_col] = 0;
                }
                for (int kk = 0; kk < not_unique.cur_pos; kk += 2) {
                    const int endpoint = not_unique[kk + 1];
                    const int endpoint_col = c.vertex_to_col[endpoint];
                    const int neighbour = not_unique[kk];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    ++not_unique_analysis[endpoint_col];
                    ++not_unique_analysis[neighbour_col];
                }
                for (int kk = 0; kk < not_unique.cur_pos; kk += 2) {
                    const int neighbour = not_unique[kk];
                    const int neighbour_col = c.vertex_to_col[neighbour];
                    const int endpoint  = not_unique[kk + 1];
                    const int endpoint_col    = c.vertex_to_col[endpoint];
                    path.reset();
                    if (!color_test.get(endpoint_col)) {
                        color_test.set(endpoint_col);
                        if (not_unique_analysis[endpoint_col] == not_unique_analysis[neighbour_col] &&
                            not_unique_analysis[endpoint_col] == c.ptn[endpoint_col] + 1) {
                            // check that path endpoints dont contain duplicates
                            bool all_unique = true;
                            color_unique.reset();
                            for (int jj = 1; jj < not_unique.cur_pos; jj += 2) {
                                const int _endpoint = not_unique[jj];
                                const int _endpoint_col = c.vertex_to_col[endpoint];
                                if (_endpoint_col == endpoint_col) {
                                    if (color_unique.get(_endpoint)) {
                                        all_unique = false;
                                        break;
                                    }
                                    color_unique.set(_endpoint);
                                }
                            }
 
                            // test_vertex connects to all vertices of endpoint_col!
                            if (all_unique && color < endpoint_col) {
                                const int path_col = c.vertex_to_col[neighbour];
                                dej_assert(c.ptn[path_col] + 1 == (not_unique_analysis[endpoint_col] * color_size));
                                dej_assert(c.ptn[endpoint_col] + 1 == not_unique_analysis[endpoint_col]);
 
 
                                int endpoint_neighbour_col = -1;
 
                                for(int jjj = 0; jjj < color_size; ++jjj ) {
                                    const int col1_vertj = c.lab[color + jjj];
 
                                    // now go through and remove paths...
                                    neighbour_list.reset();
                                    for (int jj = 0; jj < g->d[col1_vertj]; ++jj) {
                                        const int _neighbour = g->e[g->v[col1_vertj] + jj];
                                        const int _neighbour_col = c.vertex_to_col[_neighbour];
                                        //const int _neighbour_col_sz = col.ptn[_neighbour_col] + 1;
 
                                        if (_neighbour_col == path_col) {
                                            neighbour_list.push_back(_neighbour);
                                            const auto other_endpoint =
                                                    walk_to_endpoint(g, _neighbour, col1_vertj, nullptr);
                                            neighbour_to_endpoint[_neighbour] = other_endpoint.first;
                                            if(endpoint_neighbour_col == -1) {
                                                endpoint_neighbour_col = c.vertex_to_col[other_endpoint.second];
                                            }
                                        }
                                    }
 
                                    neighbour_list.sort_after_map(neighbour_to_endpoint.get_array());
 
                                    // now go through and remove paths...
                                    for (int jj = 0; jj < neighbour_list.cur_pos; ++jj) {
                                        const int _neighbour = neighbour_list[jj];
 
                                        path.reset();
                                        int _endpoint = walk_to_endpoint_collect_path(g, _neighbour,
                                                                                      col1_vertj,
                                                                                      &path);
 
                                        // mark path for deletion
                                        for (int l = 0; l < path.cur_pos; ++l) {
                                            const int del_v = path[l];
                                            dej_assert(g->d[del_v] == 2);
                                            dej_assert(!del.get(del_v));
                                            del.set(del_v);
                                            dej_assert(!path_done.get(del_v));
                                            path_done.set(del_v);
                                        }
 
                                        const int vert_orig = translate_back(col1_vertj);
 
                                        // TODO removed vertices must not contain negative values, need additional check
                                        // TODO or think about how to lift these properly
                                        recovery_strings[vert_orig].reserve(
                                                recovery_strings[vert_orig].size() +
                                                path.cur_pos);
                                        for (int l = 0; l < path.cur_pos; ++l) {
                                            dej_assert(path[l] >= 0);
                                            dej_assert(path[l] < g->v_size);
                                            const int path_v_orig = translate_back(path[l]);
                                            dej_assert(path_v_orig >= 0);
                                            dej_assert(path_v_orig < domain_size);
                                            recovery_strings[vert_orig].push_back(path_v_orig);
                                            for(size_t rsi = 0; rsi < recovery_strings[path_v_orig].size(); ++rsi) {
                                                recovery_strings[vert_orig].push_back(
                                                        recovery_strings[path_v_orig][rsi]);
                                            }
                                        }
 
 
                                        // write path into recovery_string of _endpoint
                                        const int endpoint_orig = translate_back(_endpoint);
                                        recovery_strings[endpoint_orig].reserve(
                                                recovery_strings[endpoint_orig].size() +
                                                path.cur_pos);
                                        for (int l = 0; l < path.cur_pos; ++l) {
                                            dej_assert(path[l] >= 0);
                                            dej_assert(path[l] < g->v_size);
                                            const int path_v_orig = translate_back(path[l]);
                                            dej_assert(path_v_orig >= 0);
                                            dej_assert(path_v_orig < domain_size);
                                            recovery_strings[endpoint_orig].push_back(-path_v_orig);
                                            for(size_t rsi = 0; rsi < recovery_strings[path_v_orig].size(); ++rsi) {
                                                recovery_strings[endpoint_orig].push_back(
                                                        -abs(recovery_strings[path_v_orig][rsi]));
                                            }
                                        }
 
                                        dej_assert(c.vertex_to_col[_endpoint] == endpoint_col);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
 
        // reduce vertices of degree 1 and 0, outputs corresponding automorphisms
        void red_deg10_assume_cref(dejavu::sgraph *g, int *colmap, dejavu_hook* hook) {
            if (h_deact_deg1) return;
 
            g_old_v.clear();
            g_old_v.reserve(g->v_size);
            bool has_01 = false;
            for (int i = 0; i < g->v_size; ++i) {
                const int d = g->d[i];
                has_01 = has_01 || (d == 0) || (d == 1);
                g_old_v.push_back(d);
            }
            if(!has_01) return;
 
            worklist_deg0.reset();
            worklist_deg1.reset();
 
            dejavu::markset  is_parent(g->v_size);
            dejavu::worklist pair_match(g->v_size);
            dejavu::worklist parentlist(g->v_size);
            dejavu::worklist childcount(g->v_size);
            for (int i = 0; i < g->v_size; ++i)
                childcount.push_back(0);
 
            dejavu::worklist childcount_prev(g->v_size);
            for (int i = 0; i < g->v_size; ++i)
                childcount_prev.push_back(0);
 
            dejavu::worklist_t<std::pair<int, int>> stack1(g->v_size);
            dejavu::worklist map(g->v_size);
 
            dej_assert(_automorphism_supp.cur_pos == 0);
 
            g->initialize_coloring(&c, colmap);
            for (int i = 0; i < c.domain_size;) {
                const int v = c.lab[i];
                switch (g_old_v[v]) {
                    case 0:
                        worklist_deg0.push_back(v);
                        break;
                    case 1:
                        if (c.ptn[c.vertex_to_col[v]] > 0)
                            worklist_deg1.push_back(v);
                        break;
                    default:
                        break;
                }
                i += c.ptn[i] + 1;
            }
 
            while (!worklist_deg1.empty()) {
                const int v_child = worklist_deg1.pop_back();
                if (del.get(v_child))
                    continue;
                if (g_old_v[v_child] != 1)
                    continue;
 
                const int v_child_col = c.vertex_to_col[v_child];
                const int child_col_sz = c.ptn[v_child_col] + 1;
                bool is_pairs = false;
                bool permute_parents_instead = false;
                if (child_col_sz == 1) {
                    del.set(v_child);
                    continue;
                } else {
                    parentlist.reset();
                    is_parent.reset();
                    //int last_parent_child_count = -1;
                    for (int i = v_child_col; i < v_child_col + child_col_sz; ++i) {
                        int child = c.lab[i];
 
                        // search for parent
                        const int e_pos_child = g->v[child];
                        int parent = g->e[e_pos_child];
 
                        if (is_pairs && del.get(child))
                            continue;
 
                        int search_parent = 0;
                        while (del.get(parent)) {
                            ++search_parent;
                            parent = g->e[e_pos_child + search_parent];
                        }
 
                        dej_assert(is_pairs ? c.vertex_to_col[parent] == c.vertex_to_col[child] : true);
                        if (c.vertex_to_col[parent] == c.vertex_to_col[child]) {
                            is_pairs = true;
                            del.set(child);
                            del.set(parent);
                            pair_match[child] = parent;
                            pair_match[parent] = child;
                            if (parent < child) {
                                parentlist.push_back(parent);
                            } else {
                                parentlist.push_back(child);
                            }
                            continue;
                        }
 
                        del.set(child);
 
                        // save canonical info for parent
                        edge_scratch[g->v[parent] + childcount[parent]] = child;
 
                        ++childcount[parent];
 
                        if (!is_parent.get(parent)) {
                            is_parent.set(parent);
                            childcount_prev[parent] = childcount[parent] - 1;
                            parentlist.push_back(parent);
                        }
 
                        // adjust parent degree
                        g_old_v[parent] -= 1;
                        if (g_old_v[parent] == 1 && i == v_child_col) {
                            worklist_deg1.push_back(parent);
                        } else if (g_old_v[parent] == 0 && i == v_child_col) {
                            worklist_deg0.push_back(parent);
                        }
 
                        dej_assert(g_old_v[parent] >= 0);
                    }
                }
 
                if (is_pairs) {
                    for (int j = 0; j < parentlist.cur_pos; ++j) {
                        const int first_pair_parent = parentlist[j];
                        const int pair_from = first_pair_parent;
                        const int pair_to = pair_match[pair_from];
 
                        stack1.reset();
                        map.reset();
                        map.push_back(pair_from);
                        stack1.push_back(std::pair<int, int>(g->v[pair_from], g->v[pair_from] + childcount[pair_from]));
                        while (!stack1.empty()) {
                            std::pair<int, int> from_to = stack1.pop_back();
                            int from = from_to.first;
                            const int to = from_to.second;
                            for (int f = from; f < to; ++f) {
                                const int next = edge_scratch[f];
                                //const int next = g->e[f];
                                const int from_next = g->v[next];
                                const int to_next = g->v[next] + childcount[next];
                                map.push_back(next);
                                dej_assert(next != pair_to);
                                if (from_next != to_next)
                                    stack1.push_back(std::pair<int, int>(from_next, to_next));
                            }
                        }
 
                        multiply_to_group_size(2);
 
                        dej_assert(c.vertex_to_col[pair_from] == c.vertex_to_col[pair_to]);
                        dej_assert(pair_from != pair_to);
                        int pos = 0;
 
                        //automorphism_supp.reset();
                        // descending shared_tree of child_to while writing automorphism
                        stack1.reset();
                        dej_assert(map[pos] != pair_to);
                        const int v_to_1   = pair_to;
                        const int v_from_1 = map[pos];
                        dej_assert(automorphism[v_to_1] == v_to_1);
                        dej_assert(automorphism[v_from_1] == v_from_1);
 
                        _automorphism[v_from_1] = v_to_1;
                        _automorphism[v_to_1] = v_from_1;
                        _automorphism_supp.push_back(v_from_1);
                        _automorphism_supp.push_back(v_to_1);
 
                        ++pos;
                        // child_to and child_from could have canonical strings when translated back
                        dej_assert(childcount[pair_to] == childcount[pair_from]);
                        stack1.push_back(std::pair<int, int>(g->v[pair_to], g->v[pair_to] + childcount[pair_to]));
                        while (!stack1.empty()) {
                            std::pair<int, int> from_to = stack1.pop_back();
                            int from = from_to.first;
                            const int to = from_to.second;
                            for (int f = from; f < to; ++f) {
                                const int next = edge_scratch[f];
                                //const int next = g->e[f];
                                const int from_next = g->v[next];
                                const int to_next = g->v[next] + childcount[next];
                                ++from;
                                dej_assert(next >= 0);
                                dej_assert(next < g->v_size);
                                dej_assert(map[pos] != next);
 
                                const int v_to_2 = next;
                                const int v_from_2 = map[pos];
 
                                dej_assert(_automorphism[v_to_2] == v_to_2);
                                dej_assert(_automorphism[v_from_2] == v_from_2);
                                _automorphism[v_from_2] = v_to_2;
                                _automorphism[v_to_2] = v_from_2;
                                _automorphism_supp.push_back(v_from_2);
                                _automorphism_supp.push_back(v_to_2);
 
                                ++pos;
                                if (from_next != to_next) // there was a semicolon here, should have been bug
                                    stack1.push_back(std::pair<int, int>(from_next, to_next));
                            }
                        }
 
                        dej_assert(pos == map.cur_pos);
 
                        dej_assert(g_old_v[pair_to] == 1);
                        dej_assert(g_old_v[pair_from] == 1);
                        dej_assert(del.get(pair_to));
                        dej_assert(del.get(pair_from));
 
                        pre_hook(g->v_size, _automorphism.get_array(), _automorphism_supp.cur_pos,
                                 _automorphism_supp.get_array(), hook);
 
                        reset_automorphism(_automorphism.get_array(), _automorphism_supp.cur_pos,
                                           _automorphism_supp.get_array());
                        _automorphism_supp.reset();
                    }
 
                    for (int j = 0; j < parentlist.cur_pos; ++j) {
                        const int first_pair_parent = parentlist[j];
                        const int pair_from = first_pair_parent;
                        const int pair_to = pair_match[pair_from];
 
                        const int original_parent = translate_back(first_pair_parent);
                        const int original_pair_to = translate_back(pair_to);
                        recovery_strings[original_parent].push_back(original_pair_to);
                        for (size_t s = 0; s < recovery_strings[original_pair_to].size(); ++s)
                            recovery_strings[original_parent].push_back(
                                    recovery_strings[original_pair_to][s]);
 
                        // also write stack of original_pair_to to canonical recovery string of original_parent, since
                        // recovery_strings has not been updated with progress made in this routine
                        stack1.reset();
                        stack1.push_back(std::pair<int, int>(g->v[pair_to], g->v[pair_to] + childcount[pair_to]));
                        while (!stack1.empty()) {
                            std::pair<int, int> from_to = stack1.pop_back();
                            int from = from_to.first;
                            const int to = from_to.second;
                            for (int f = from; f < to; ++f) {
                                const int next = edge_scratch[f];
                                //const int next = g->e[f];
                                const int from_next = g->v[next];
                                const int to_next = g->v[next] + childcount[next];
                                const int orig_next = translate_back(next);
                                recovery_strings[original_parent].push_back(orig_next);
                                for (size_t s = 0; s < recovery_strings[orig_next].size(); ++s)
                                    recovery_strings[original_parent].push_back(
                                            recovery_strings[orig_next][s]);
                                if (from_next != to_next)
                                    stack1.push_back(std::pair<int, int>(from_next, to_next));
                            }
                        }
                    }
 
                    permute_parents_instead = true;
                }
 
                while (!parentlist.empty()) {
                    int parent, childcount_from, childcount_to, child_from;
                    if (!permute_parents_instead) {
                        parent = parentlist.pop_back();
                        childcount_from = childcount_prev[parent];
                        childcount_to = childcount[parent];
                    } else {
                        parent = -1;
                        childcount_from = 0;
                        childcount_to = parentlist.cur_pos;
                    }
                    // automorphism 1: long cycle (c1 ... cn)
                    dej_assert(childcount_to - childcount_from > 0);
                    if (childcount_to - childcount_from == 1) {
                        if (permute_parents_instead)
                            break;
                        continue;
                    }
                    if (!permute_parents_instead) {
                        child_from = edge_scratch[g->v[parent] + childcount_from];
                        //child_from = g->e[g->v[parent] + childcount_from];
                    } else {
                        child_from = parentlist[0];
                    }
 
                    // descending shared_tree of child_from while writing map
                    stack1.reset();
                    map.reset();
                    map.push_back(child_from);
                    stack1.push_back(std::pair<int, int>(g->v[child_from], g->v[child_from] + childcount[child_from]));
                    while (!stack1.empty()) {
                        std::pair<int, int> from_to = stack1.pop_back();
                        int from = from_to.first;
                        const int to = from_to.second;
                        for (int f = from; f < to; ++f) {
                            const int next = edge_scratch[f];
                            //const int next = g->e[f];
                            const int from_next = g->v[next];
                            const int to_next = g->v[next] + childcount[next];
                            map.push_back(next);
                            dej_assert(next != parent);
                            if (from_next != to_next)
                                stack1.push_back(std::pair<int, int>(from_next, to_next));
                        }
                    }
                    int j = 2;
                    for (int i = childcount_from + 1; i < childcount_to; ++i) {
                        multiply_to_group_size(j);
                        ++j;
                        int child_to;
                        if (!permute_parents_instead) {
                            child_to = edge_scratch[g->v[parent] + i];
                            //child_to = g->e[g->v[parent] + i];
                        } else {
                            child_to = parentlist[i];
                        }
                        dej_assert(c.vertex_to_col[child_from] == c.vertex_to_col[child_to]);
                        dej_assert(child_from != child_to);
                        int pos = 0;
 
                        _automorphism_supp.reset();
 
                        // descending shared_tree of child_to while writing automorphism
                        stack1.reset();
                        dej_assert(map[pos] != child_to);
 
                        const int v_to_1 = child_to;
                        const int v_from_1 = map[pos];
                        dej_assert(_automorphism[v_to_1] == v_to_1);
                        dej_assert(_automorphism[v_from_1] == v_from_1);
                        _automorphism[v_from_1] = v_to_1;
                        _automorphism[v_to_1] = v_from_1;
                        _automorphism_supp.push_back(v_from_1);
                        _automorphism_supp.push_back(v_to_1);
 
                        ++pos;
                        // child_to and child_from could have canonical strings when translated back
                        dej_assert(childcount[child_to] == childcount[child_from]);
                        stack1.push_back(std::pair<int, int>(g->v[child_to], g->v[child_to] + childcount[child_to]));
                        while (!stack1.empty()) {
                            std::pair<int, int> from_to = stack1.pop_back();
                            int from = from_to.first;
                            const int to = from_to.second;
                            for (int f = from; f < to; ++f) {
                                const int next = edge_scratch[f];
                                //const int next      = g->e[f];
                                const int from_next = g->v[next];
                                const int to_next = g->v[next] + childcount[next];
                                ++from;
                                dej_assert(next >= 0);
                                dej_assert(next < g->v_size);
                                dej_assert(map[pos] != next);
 
                                const int v_to_2 = next;
                                const int v_from_2 = map[pos];
                                dej_assert(_automorphism[v_to_2] == v_to_2);
                                dej_assert(_automorphism[v_from_2] == v_from_2);
                                _automorphism[v_from_2] = v_to_2;
                                _automorphism[v_to_2] = v_from_2;
                                _automorphism_supp.push_back(v_from_2);
                                _automorphism_supp.push_back(v_to_2);
 
                                ++pos;
                                if (from_next != to_next) // there was a semicolon here, should have been bug
                                    stack1.push_back(std::pair<int, int>(from_next, to_next));
                            }
                        }
 
                        dej_assert(pos == map.cur_pos);
                        dej_assert(g_old_v[child_to] == 1);
                        dej_assert(g_old_v[child_from] == 1);
                        dej_assert(del.get(child_to));
                        dej_assert(del.get(child_from));
 
                        pre_hook(g->v_size, _automorphism.get_array(), _automorphism_supp.cur_pos,
                                 _automorphism_supp.get_array(), hook);
 
                        reset_automorphism(_automorphism.get_array(), _automorphism_supp.cur_pos,
                                           _automorphism_supp.get_array());
                        _automorphism_supp.reset();
                    }
                    if (permute_parents_instead) {
                        break;      // :-)
                    }
                }
                parentlist.reset();
            }
 
            if (g->v_size > 1) {
                // search for parents that still remain in the graph, and rewrite childlist structure into canonical string
                for (int _i = 0; _i < g->v_size; ++_i) {
                    if (!del.get(_i) && childcount[_i] > 0 && c.ptn[c.vertex_to_col[_i]] > 0) {
                        // should it be childcount[_i] > 0 or childcount[_i] >= 0? was childcount[_i] >= 0 but that doesnt make sense?
                        stack1.reset();
                        stack1.push_back(std::pair<int, int>(g->v[_i], g->v[_i] + childcount[_i]));
                        const int orig_i = translate_back(_i);
                        recovery_strings[orig_i].reserve(
                                recovery_strings[orig_i].size() + childcount[_i]);
                        while (!stack1.empty()) {
                            std::pair<int, int> from_to = stack1.pop_back();
                            int from = from_to.first;
                            const int to = from_to.second;
                            if (from == to) {
                                continue;
                            } else {
                                const int next = edge_scratch[from];
                                //const int next = g->e[from];
                                const int from_next = g->v[next];
                                const int to_next = g->v[next] + childcount[next];
                                ++from;
                                const int orig_next = translate_back(next);
                                recovery_strings[orig_i].push_back(orig_next);
                                // write canonical recovery string of orig_next into orig_i, since it is now represented by
                                // orig_i
                                dej_assert(next != _i);
                                dej_assert(orig_next != orig_i);
                                for (size_t j = 0; j < recovery_strings[orig_next].size(); ++j)
                                    recovery_strings[orig_i].push_back(
                                            recovery_strings[orig_next][j]);
                                stack1.push_back(std::pair<int, int>(from, to));
                                stack1.push_back(std::pair<int, int>(from_next, to_next));
                            }
                        }
                    }
                }
            }
 
 
            while (!worklist_deg0.empty()) {
                const int v_child = worklist_deg0.pop_back();
                if (del.get(v_child))
                    continue;
                dej_assert(g_old_v[v_child] == 0);
 
                is_parent.reset();
                const int v_child_col = c.vertex_to_col[v_child];
                const int child_col_sz = c.ptn[v_child_col] + 1;
                const int parent_from = c.lab[v_child_col];
                del.set(parent_from);
 
                if (child_col_sz == 1)
                    continue;
                int j = 2;
                for (int i = v_child_col + 1; i < v_child_col + child_col_sz; ++i) {
                    multiply_to_group_size(j);
                    ++j;
                    const int parent_to = c.lab[i];
                    dej_assert(g_old_v[parent_to] == 0);
                    del.set(parent_to);
 
                    dej_assert(parent_from != parent_to);
 
                    dej_assert(_automorphism_supp.cur_pos == 0);
 
                    _automorphism[parent_to] = parent_from;
                    _automorphism[parent_from] = parent_to;
                    _automorphism_supp.push_back(parent_from);
                    _automorphism_supp.push_back(parent_to);
 
                    pre_hook(g->v_size, _automorphism.get_array(), _automorphism_supp.cur_pos,
                             _automorphism_supp.get_array(), hook);
 
                    reset_automorphism(_automorphism.get_array(), _automorphism_supp.cur_pos,
                                       _automorphism_supp.get_array());
                    _automorphism_supp.reset();
                }
            }
        }
 
        // perform edge flips according to quotient graph
        void red_quotient_edge_flip(dejavu::sgraph *g, int *colmap) { // TODO could still optimize further ...
            if (g->v_size <= 1) return;
 
            del_e.reset();
            worklist_deg0.reset(); // serves as the "test vector"
 
            //coloring test_col;
            g->initialize_coloring(&c, colmap);
 
            for (int y = 0; y < g->v_size; ++y) worklist_deg0[y] = 0;
 
            for (int i = 0; i < g->v_size;) {
                int v = c.lab[i];
                for (int f = g->v[v]; f < g->v[v] + g->d[v]; ++f) {
                    const int v_neigh = g->e[f];
                    worklist_deg0[c.vertex_to_col[v_neigh]] += 1;
                }
 
                for (int ii = 0; ii < c.ptn[i] + 1; ++ii) {
                    const int vx = c.lab[i + ii];
                    bool skipped_none = true;
                    for (int f = g->v[vx]; f < g->v[vx] + g->d[vx]; ++f) {
                        const int v_neigh = g->e[f];
                        if (worklist_deg0[c.vertex_to_col[v_neigh]] ==
                                c.ptn[c.vertex_to_col[v_neigh]] + 1) {
                            del_e.set(f); // mark edge for deletion (reverse edge is handled later automatically)
                            skipped_none = false;
                        }
                    }
                    if(skipped_none) break;
                }
 
                for (int f = g->v[v]; f < g->v[v] + g->d[v]; ++f) {
                    const int v_neigh = g->e[f];
                    worklist_deg0[c.vertex_to_col[v_neigh]] = 0;
                }
 
                i += c.ptn[i] + 1;
            }
        }
 
        static void copy_coloring_to_colmap(const coloring *c, int *colmap) {
            for (int i = 0; i < c->domain_size; ++i) {
                colmap[i] = c->vertex_to_col[i];
            }
        }
 
        // deletes vertices marked in 'del'
        // assumes that g->v points to g->e in ascending order
        void perform_del(dejavu::sgraph *g, int *colmap) {
            // copy some stuff
            g_old_v.clear();
            translate_layer_fwd.clear();
            translate_layer_bwd.clear();
 
            translate_layer_bwd.reserve(backward_translation_layers[backward_translation_layers.size() - 1].size());
            translate_layer_fwd.reserve(g->v_size);
            for (size_t i = 0; i < backward_translation_layers[backward_translation_layers.size() - 1].size(); ++i)
                translate_layer_bwd.push_back(
                        backward_translation_layers[backward_translation_layers.size() - 1][i]); // TODO this is shit
 
            // create translation array from old graph to new graph vertices
            int cnt = 0;
            int new_vsize = 0;
            // int pos_now = 0;
            for (int i = 0; i < g->v_size; ++i) {
                if (!del.get(i)) {
                    translate_layer_fwd.push_back(cnt);
                    const int translate_back = translate_layer_bwd[translate_layer_fwd.size() - 1];
                    backward_translation_layers[backward_translation_layers.size() - 1][cnt] = translate_back;
                    ++cnt;
                    ++new_vsize;
                } else {
                    translate_layer_fwd.push_back(-1);
                }
            }
 
            if (new_vsize == 0 || new_vsize == 1) {
                g->v_size = 0;
                g->e_size = 0;
                return;
            }
 
            backward_translation_layers[backward_translation_layers.size() - 1].resize(cnt);
 
            g_old_v.reserve(g->v_size);
            for (int i = 0; i < g->v_size; ++i) {
                g_old_v.push_back(g->v[i]);
            }
 
            // make graph smaller using the translation array
            int epos = 0;
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
 
                if (new_v >= 0) {
                    int new_d = 0;
                    dej_assert(new_v < new_vsize);
                    g->v[new_v] = epos;
                    for (int j = g_old_v[old_v]; j < g_old_v[old_v] + g->d[old_v]; ++j) {
                        const int ve = g->e[j];                          // assumes ascending order!
                        const int new_ve = translate_layer_fwd[ve];
                        if (new_ve >= 0) {
                            dej_assert(new_ve < new_vsize);
                            dej_assert(new_ve >= 0);
                            ++new_d;
                            dej_assert(j >= epos);
                            g->e[epos] = new_ve;                         // assumes ascending order!
                            ++epos;
                        }
                    }
 
                    g_old_v[old_v] = new_d;
                }
            }
 
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
                dej_assert(old_v >= new_v);
                if (new_v >= 0) {
                    dej_assert(old_v >= 0);
                    dej_assert(new_v >= 0);
                    dej_assert(old_v < g->v_size);
                    dej_assert(new_v < new_vsize);
                    colmap[new_v] = colmap[old_v];
                }
            }
 
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
                if (new_v >= 0) {
                    g->d[new_v] = g_old_v[old_v];
                }
            }
 
            dej_assert(new_vsize <= g->v_size);
            dej_assert(epos <= g->e_size);
 
            g->e_size = epos;
            g->v_size = new_vsize;
            del.reset();
        }
 
        // deletes vertices marked in 'del'
        // assumes that g->v points to g->e in ascending order
        void perform_del_edge(dejavu::sgraph *g) {
            if (g->v_size <= 1)
                return;
 
            // int pre_esize = g->e_size;
            // copy some stuff
            g_old_v.clear();
            g_old_v.reserve(g->v_size);
 
            for (int i = 0; i < g->v_size; ++i) {
                g_old_v.push_back(g->v[i]);
            }
 
            // create translation array from old graph to new graph vertices
            // int cnt = 0;
            int new_vsize = g->v_size;
 
            // make graph smaller using the translation array
            int epos = 0;
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = old_v;
 
                if (new_v >= 0) {
                    int new_d = 0;
                    g->v[new_v] = epos;
                    for (int j = g_old_v[old_v]; j < g_old_v[old_v] + g->d[old_v]; ++j) {
                        const int ve = g->e[j];
                        const int new_ve = ve;
                        if (!del_e.get(j)) {
                            dej_assert(new_ve < new_vsize);
                            dej_assert(new_ve >= 0);
                            ++new_d;
                            g->e[epos] = new_ve;
                            ++epos;
                        }
                    }
 
                    g->d[new_v] = new_d;
                }
            }
 
            g->e_size = epos;
            g->v_size = new_vsize;
        }
 
        // deletes vertices marked in 'del'
        // for vertices v where add_edge_buff_act[v] is set, in turn adds edges add_edge_buff_act[v]
        // assumes that g->v points to g->e in ascending order
        // assumes that degree of a vertex stays the same or gets smaller
        void perform_del_add_edge(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1)
                return;
 
            // copy some stuff
            g_old_v.clear();
            translate_layer_fwd.clear();
            translate_layer_bwd.clear();
 
            for (size_t i = 0; i < backward_translation_layers[backward_translation_layers.size() - 1].size(); ++i)
                translate_layer_bwd.push_back(backward_translation_layers[backward_translation_layers.size() - 1][i]);
 
            // create translation array from old graph to new graph vertices
            int cnt = 0;
            int new_vsize = 0;
            for (int i = 0; i < g->v_size; ++i) {
                worklist_deg0[i] = colmap[i];
                if (!del.get(i)) {
                    translate_layer_fwd.push_back(cnt);
                    const int translate_back = translate_layer_bwd[translate_layer_fwd.size() - 1];
                    backward_translation_layers[backward_translation_layers.size() - 1][cnt] = translate_back;
                    ++cnt;
                    ++new_vsize;
                } else {
                    //translation_layers[fwd_ind].push_back(-1);
                    translate_layer_fwd.push_back(-1);
                }
            }
 
            if (new_vsize == g->v_size)
                return;
 
            if (new_vsize == 0 || new_vsize == 1) {
                g->v_size = 0;
                g->e_size = 0;
                return;
            }
 
            g_old_v.reserve(g->v_size);
 
            for (int i = 0; i < g->v_size; ++i) {
                g_old_v.push_back(g->v[i]);
            }
 
            backward_translation_layers[backward_translation_layers.size() - 1].resize(cnt);
 
            // make graph smaller using the translation array
            int epos = 0;
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[old_v];
 
                if (new_v >= 0) {
                    int new_d = 0;
                    g->v[new_v] = epos;
                    for (int j = g_old_v[old_v]; j < g_old_v[old_v] + g->d[old_v]; ++j) {
                        const int ve = g->e[j];
                        const int new_ve = translate_layer_fwd[ve];
                        if (new_ve >= 0) {
                            dej_assert(new_ve < new_vsize);
                            dej_assert(new_ve >= 0);
                            ++new_d;
                            g->e[epos] = new_ve;
                            ++epos;
                        }
                    }
                    if (add_edge_buff_act.get(old_v)) {
                        while (!add_edge_buff[old_v].empty()) {
                            const int edge_to_old = add_edge_buff[old_v].back();
                            add_edge_buff[old_v].pop_back();
                            //const int edge_to_old = add_edge_buff[old_v];
                            dej_assert(add_edge_buff_act.get(edge_to_old));
                            //const int edge_to_new = translation_layers[fwd_ind][edge_to_old];
                            const int edge_to_new = translate_layer_fwd[edge_to_old];
                            dej_assert(edge_to_old >= 0);
                            dej_assert(edge_to_old <= domain_size);
                            dej_assert(edge_to_new >= 0);
                            dej_assert(edge_to_new <= new_vsize);
                            ++new_d;
                            g->e[epos] = edge_to_new;
                            ++epos;
                        }
                    }
                    g_old_v[old_v] = new_d;
                }
            }
 
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
                if (new_v >= 0) {
                    g->d[new_v] = g_old_v[old_v];
                }
            }
 
            g->e_size = epos;
 
            // adapt colmap for remaining vertices
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                //const int new_v = translation_layers[fwd_ind][i];
                const int new_v = translate_layer_fwd[i];
                dej_assert(new_v < new_vsize);
                if (new_v >= 0) {
                    dej_assert(colmap[old_v] >= 0);
                    //assert(colmap[old_v] < domain_size);
                    //colmap[new_v] = colmap[old_v];
                    colmap[new_v] = worklist_deg0[old_v];
                }
            }
 
            g->v_size = new_vsize;
 
            for (int i = 0; i < g->v_size; ++i) {
                dej_assert(g->d[i] > 0 ? g->v[i] < g->e_size : true);
                dej_assert(g->d[i] > 0 ? g->v[i] >= 0 : true);
                dej_assert(g->d[i] >= 0);
                dej_assert(g->d[i] < g->v_size);
            }
            for (int i = 0; i < g->e_size; ++i) {
                dej_assert(g->e[i] < g->v_size);
                dej_assert(g->e[i] >= 0);
            }
 
            add_edge_buff_act.reset();
            del.reset();
        }
 
        // deletes vertices marked in 'del'
        // for vertices v where add_edge_buff_act[v] is set, in turn adds edges add_edge_buff_act[v]
        // assumes that g->v points to g->e in ascending order
        // assumes that degree of a vertex stays the same or gets smaller
        void perform_del_add_edge_general(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1)
                return;
 
            // copy some stuff
            g_old_v.clear();
            g_old_e.clear();
            translate_layer_fwd.clear();
            translate_layer_bwd.clear();
 
            for (size_t i = 0; i < backward_translation_layers[backward_translation_layers.size() - 1].size(); ++i)
                translate_layer_bwd.push_back(backward_translation_layers[backward_translation_layers.size() - 1][i]);
 
            // create translation array from old graph to new graph vertices
            int cnt = 0;
            int new_vsize = 0;
            for (int i = 0; i < g->v_size; ++i) {
                worklist_deg0[i] = colmap[i];
                if (!del.get(i)) {
                    translate_layer_fwd.push_back(cnt);
                    const int translate_back = translate_layer_bwd[translate_layer_fwd.size() - 1];
                    backward_translation_layers[backward_translation_layers.size() - 1][cnt] = translate_back;
                    ++cnt;
                    ++new_vsize;
                } else {
                    //translation_layers[fwd_ind].push_back(-1);
                    translate_layer_fwd.push_back(-1);
                }
            }
 
            if (new_vsize == g->v_size)
                return;
 
            if (new_vsize == 0 || new_vsize == 1) {
                g->v_size = 0;
                g->e_size = 0;
                return;
            }
 
            g_old_v.reserve(g->v_size);
 
            for (int i = 0; i < g->v_size; ++i) {
                g_old_v.push_back(g->v[i]);
            }
            for (int i = 0; i < g->e_size; ++i) {
                g_old_e.push_back(g->e[i]);
            }
 
            backward_translation_layers[backward_translation_layers.size() - 1].resize(cnt);
 
            // make graph smaller using the translation array
            int epos = 0;
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[old_v];
 
                if (new_v >= 0) {
                    int new_d = 0;
                    g->v[new_v] = epos;
                    for (int j = g_old_v[old_v]; j < g_old_v[old_v] + g->d[old_v]; ++j) {
                        const int ve = g_old_e[j];
                        const int new_ve = translate_layer_fwd[ve];
                        if (new_ve >= 0) {
                            dej_assert(new_ve < new_vsize);
                            dej_assert(new_ve >= 0);
                            ++new_d;
                            g->e[epos] = new_ve;
                            ++epos;
                        }
                    }
                    if (add_edge_buff_act.get(old_v)) {
                        while (!add_edge_buff[old_v].empty()) {
                            const int edge_to_old = add_edge_buff[old_v].back();
                            add_edge_buff[old_v].pop_back();
                            //const int edge_to_old = add_edge_buff[old_v];
                            dej_assert(add_edge_buff_act.get(edge_to_old));
                            //const int edge_to_new = translation_layers[fwd_ind][edge_to_old];
                            const int edge_to_new = translate_layer_fwd[edge_to_old];
                            dej_assert(edge_to_old >= 0);
                            dej_assert(edge_to_old <= domain_size);
                            dej_assert(edge_to_new >= 0);
                            dej_assert(edge_to_new <= new_vsize);
                            ++new_d;
                            g->e[epos] = edge_to_new;
                            ++epos;
                        }
                    }
                    g_old_v[old_v] = new_d;
                }
            }
 
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
                if (new_v >= 0) {
                    g->d[new_v] = g_old_v[old_v];
                }
            }
 
            g->e_size = epos;
 
            // adapt colmap for remaining vertices
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                //const int new_v = translation_layers[fwd_ind][i];
                const int new_v = translate_layer_fwd[i];
                dej_assert(new_v < new_vsize);
                if (new_v >= 0) {
                    dej_assert(colmap[old_v] >= 0);
                    //assert(colmap[old_v] < domain_size);
                    //colmap[new_v] = colmap[old_v];
                    colmap[new_v] = worklist_deg0[old_v];
                }
            }
 
            g->v_size = new_vsize;
 
            for (int i = 0; i < g->v_size; ++i) {
                dej_assert(g->d[i] > 0 ? g->v[i] < g->e_size : true);
                dej_assert(g->d[i] > 0 ? g->v[i] >= 0 : true);
                dej_assert(g->d[i] >= 0);
                dej_assert(g->d[i] < g->v_size);
            }
            for (int i = 0; i < g->e_size; ++i) {
                dej_assert(g->e[i] < g->v_size);
                dej_assert(g->e[i] >= 0);
            }
 
            add_edge_buff_act.reset();
            del.reset();
        }
 
        // marks all discrete vertices in 'del'
        void mark_discrete_for_deletion(dejavu::sgraph *g, int *colmap) {
            //int discrete_cnt = 0;
            worklist_deg0.reset();
            for (int i = 0; i < domain_size; ++i) {
                worklist_deg0.push_back(0);
            }
            for (int i = 0; i < g->v_size; ++i) {
                dej_assert(colmap[i] < domain_size);
                worklist_deg0[colmap[i]]++;
            }
            for (int i = 0; i < g->v_size; ++i) {
                if (worklist_deg0[colmap[i]] == 1)
                    del.set(i);
            }
            worklist_deg0.reset();
        }
 
        // deletes all discrete vertices
        void perform_del_discrete(dejavu::sgraph *g, int *colmap) {
            if (g->v_size <= 1)
                return;
 
            int discrete_cnt = 0;
            worklist_deg0.reset();
            for (int i = 0; i < domain_size; ++i) {
                worklist_deg0.push_back(0);
            }
            for (int i = 0; i < g->v_size; ++i) {
                dej_assert(colmap[i] < domain_size);
                worklist_deg0[colmap[i]]++;
            }
            for (int i = 0; i < g->v_size; ++i) {
                discrete_cnt += (worklist_deg0[colmap[i]] == 1);
            }
            if (discrete_cnt == g->v_size) {
                g->v_size = 0;
                g->e_size = 0;
                return;
            }
            if (discrete_cnt == 0) {
                return;
            }
 
            // copy some stuff
            g_old_v.clear();
            translate_layer_fwd.clear();
            translate_layer_bwd.clear();
 
            for (size_t i = 0; i < backward_translation_layers[backward_translation_layers.size() - 1].size(); ++i)
                translate_layer_bwd.push_back(backward_translation_layers[backward_translation_layers.size() - 1][i]);
 
            g_old_v.reserve(g->v_size);
 
            for (int i = 0; i < g->v_size; ++i) {
                g_old_v.push_back(g->v[i]);
            }
 
            // create translation array from old graph to new graph vertices
            int cnt = 0;
            int new_vsize = 0;
            for (int i = 0; i < g->v_size; ++i) {
                if (worklist_deg0[colmap[i]] != 1) {
                    translate_layer_fwd.push_back(cnt);
                    const int translate_back = translate_layer_bwd[translate_layer_fwd.size() - 1];
                    backward_translation_layers[backward_translation_layers.size() - 1][cnt] = translate_back;
                    ++cnt;
                    ++new_vsize;
                } else {
                    translate_layer_fwd.push_back(-1);
                }
            }
 
            backward_translation_layers[backward_translation_layers.size() - 1].resize(cnt);
 
            // make graph smaller using the translation array
            int epos = 0;
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
 
                if (new_v >= 0) {
                    int new_d = 0;
                    dej_assert(new_v < new_vsize);
                    g->v[new_v] = epos;
                    for (int j = g_old_v[old_v]; j < g_old_v[old_v] + g->d[old_v]; ++j) {
                        const int ve = g->e[j];                          // assumes ascending order!
                        const int new_ve = ve>=0?translate_layer_fwd[ve]:-1;
                        if (new_ve >= 0) {
                            dej_assert(new_ve < new_vsize);
                            dej_assert(new_ve >= 0);
                            ++new_d;
                            dej_assert(j >= epos);
                            g->e[epos] = new_ve;                         // assumes ascending order!
                            ++epos;
                        }
                    }
 
                    g_old_v[old_v] = new_d;
                }
            }
 
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
                dej_assert(old_v >= new_v);
                if (new_v >= 0) {
                    dej_assert(old_v >= 0);
                    dej_assert(new_v >= 0);
                    dej_assert(old_v < g->v_size);
                    dej_assert(new_v < new_vsize);
                    colmap[new_v] = colmap[old_v];
                }
            }
 
            for (int i = 0; i < g->v_size; ++i) {
                const int old_v = i;
                const int new_v = translate_layer_fwd[i];
                if (new_v >= 0) {
                    g->d[new_v] = g_old_v[old_v];
                }
            }
 
            dej_assert(new_vsize <= g->v_size);
            dej_assert(epos <= g->e_size);
 
            g->e_size = epos;
            g->v_size = new_vsize;
            del.reset();
        }
 
        // given automorphism of reduced graph, reconstructs automorphism of the original graph
        // does not optimize for consecutive calls
        void pre_hook(int, const int *_aut, int _supp, const int *_aut_supp, dejavu_hook* hook) {
            if(hook == nullptr)
                return;
 
            automorphism_supp.reset();
 
            bool use_aux_auto = false;
 
            for (int i = 0; i < _supp; ++i) {
                const int v_from = _aut_supp[i];
                const int orig_v_from = translate_back(v_from);
                const int v_to = _aut[v_from];
                dej_assert(v_from != v_to);
                const int orig_v_to = translate_back(v_to);
                dej_assert(v_from >= 0);
                dej_assert(v_to >= 0);
                dej_assert(orig_v_from < domain_size);
                dej_assert(orig_v_from >= 0);
                dej_assert(orig_v_to < domain_size);
                dej_assert(orig_v_to >= 0);
 
                const bool added_vertex = !recovery_strings[orig_v_from].empty() &&
                                           recovery_strings[orig_v_from][0] == INT32_MAX;
 
                dej_assert(automorphism[orig_v_from] == orig_v_from || added_vertex);
 
                if(!added_vertex) {
                    automorphism[orig_v_from] = orig_v_to;
                    automorphism_supp.push_back(orig_v_from);
                } else {
                    dej_assert(recovery_strings[orig_v_to].size() == 1);
                }
 
                dej_assert(recovery_strings[orig_v_to].size() == recovery_strings[orig_v_from].size());
 
                for (size_t j = 0; j < recovery_strings[orig_v_to].size(); ++j) {
                    const int v_from_t = recovery_strings[orig_v_from][j];
                    const int v_to_t   = recovery_strings[orig_v_to][j];
                    dej_assert((v_from_t == INT32_MAX) == (v_to_t == INT32_MAX));
                    if(v_from_t == INT32_MAX) continue;
 
                    dej_assert(v_from_t >  0?v_to_t >= 0:true);
                    dej_assert(v_to_t   >  0?v_from_t >= 0:true);
                    dej_assert(v_from_t <  0?v_to_t <= 0:true);
                    dej_assert(v_to_t   <  0?v_from_t <= 0:true);
                    if(v_from_t >= 0 && v_to_t >= 0) {
                        dej_assert(automorphism[v_from_t] == v_from_t);
                        automorphism[v_from_t] = v_to_t;
                        automorphism_supp.push_back(v_from_t);
                    } else {
                        const int abs_v_from_t = abs(v_from_t);
                        const int abs_v_to_t   = abs(v_to_t);
 
                        dej_assert(aux_automorphism[abs_v_from_t] == abs_v_from_t);
                        aux_automorphism[abs_v_from_t] = abs_v_to_t;
                        aux_automorphism_supp.push_back(abs_v_from_t);
 
                        use_aux_auto = true;
                    }
                }
            }
 
            if(use_aux_auto) {
                for(int i = 0; i < aux_automorphism_supp.cur_pos; ++i) {
                    const int v_from = aux_automorphism_supp[i];
                    before_move[v_from]  = automorphism[aux_automorphism[v_from]];
                }
                for(int i = 0; i < aux_automorphism_supp.cur_pos; ++i) {
                    const int v_from = aux_automorphism_supp[i];
                    if(automorphism[v_from] == v_from) {
                        automorphism_supp.push_back(v_from);
                    }
                    dej_assert(automorphism[v_from] != before_move[v_from]);
                    automorphism[v_from] = before_move[v_from];
                }
                reset_automorphism(aux_automorphism.get_array(), aux_automorphism_supp.cur_pos,
                                   aux_automorphism_supp.get_array());
                aux_automorphism_supp.reset();
            }
 
            const int end_pos = automorphism_supp.cur_pos;
 
            if(!recovery_edge_attached.empty()) {
                for (int i = 0; i < end_pos; ++i) {
                    recovery_attached_edges_to_automorphism(automorphism_supp[i], automorphism.get_array(),
                                                            &automorphism_supp, &aux_automorphism, &before_move);
                }
            }
 
            if(hook != nullptr)
                (*hook)(domain_size, automorphism.get_array(), automorphism_supp.cur_pos, automorphism_supp.get_array());
            reset_automorphism(automorphism.get_array(), automorphism_supp.cur_pos, automorphism_supp.get_array());
            automorphism_supp.reset();
        }
 
    public:
        // given automorphism of reduced graph, reconstructs automorphism of the original graph
        void
        pre_hook_buffered(int _n, const int *_aut, int _supp, const int *_aut_supp, dejavu_hook* hook) {
            if(hook == nullptr) {
                return;
            }
 
            meld_translation_layers();
            // TODO what this should really do is bake recovery_strings into a single string, and edge attachmenets to
            // TODO a proper non-linked graph structure
 
            automorphism_supp.reset();
            bool use_aux_auto = false;
 
            if(_supp >= 0) {
                for (int i = 0; i < _supp; ++i) {
                    const int _v_from = _aut_supp[i];
                    const int v_from  = decomposer==nullptr?_v_from:decomposer->map_back(current_component, _v_from);
                    dej_assert(v_from >= 0 && v_from < domain_size);
                    const int orig_v_from = backward_translation[v_from];
                    const int _v_to = _aut[_v_from];
                    const int v_to  = decomposer==nullptr?_v_to:decomposer->map_back(current_component, _v_to);
                    dej_assert(v_from != v_to);
                    const int orig_v_to = backward_translation[v_to];
                    dej_assert((unsigned int)v_from < backward_translation.size());
                    dej_assert(v_from >= 0);
                    dej_assert((unsigned int)v_to < backward_translation.size());
                    dej_assert(v_to >= 0);
                    dej_assert(orig_v_from < domain_size);
                    dej_assert(orig_v_from >= 0);
                    dej_assert(orig_v_to < domain_size);
                    dej_assert(orig_v_to >= 0);
                    dej_assert(automorphism[orig_v_from] == orig_v_from);
                    const int to_recovery_string_start_pt = baked_recovery_string_pt[v_to].first;
                    const int to_recovery_string_end_pt   = baked_recovery_string_pt[v_to].second;
                    const int to_recovery_string_size     = to_recovery_string_end_pt - to_recovery_string_start_pt;
                    const int from_recovery_string_start_pt = baked_recovery_string_pt[v_from].first;
                    if(to_recovery_string_size == 0 || baked_recovery_string[to_recovery_string_start_pt] != INT32_MAX) {
                        automorphism[orig_v_from] = orig_v_to;
                        automorphism_supp.push_back(orig_v_from);
                    }
 
                    dej_assert((baked_recovery_string_pt[v_from].second - from_recovery_string_start_pt) ==
                               to_recovery_string_size);
 
                    for (int j = 0; j < to_recovery_string_size; ++j) {
                        const int v_from_t = baked_recovery_string[from_recovery_string_start_pt + j];
                        const int v_to_t   = baked_recovery_string[to_recovery_string_start_pt + j];
 
                        if(v_from_t == INT32_MAX) continue;
 
                        dej_assert(abs(v_from_t) >= 0);
                        dej_assert(abs(v_from_t) <  domain_size);
                        dej_assert(abs(v_to_t) >= 0);
                        dej_assert(abs(v_to_t) <  domain_size);
 
                        dej_assert(v_from_t >  0?v_to_t >= 0:true);
                        dej_assert(v_to_t   >  0?v_from_t >= 0:true);
                        dej_assert(v_from_t <  0?v_to_t <= 0:true);
                        dej_assert(v_to_t   <  0?v_from_t <= 0:true);
                        if(v_from_t >= 0 && v_to_t >= 0) {
                            dej_assert(automorphism[v_from_t] == v_from_t);
                            automorphism[v_from_t] = v_to_t;
                            automorphism_supp.push_back(v_from_t);
                        } else {
                            const int abs_v_from_t = abs(v_from_t);
                            const int abs_v_to_t   = abs(v_to_t);
 
                            aux_automorphism[abs_v_from_t] = abs_v_to_t;
                            aux_automorphism_supp.push_back(abs_v_from_t);
 
                            use_aux_auto = true;
                        }
                    }
                }
            } else {
                for (int i = 0; i < _n; ++i) {
                    const int _v_from = i;
                    const int v_from  = decomposer==nullptr?_v_from:decomposer->map_back(current_component, _v_from);
                    const int orig_v_from = backward_translation[v_from];
                    const int _v_to = _aut[_v_from];
                    const int v_to  = decomposer==nullptr?_v_to:decomposer->map_back(current_component, _v_to);
                    if(v_from == v_to)
                        continue;
                    const int orig_v_to = backward_translation[v_to];
                    dej_assert((unsigned int)v_from < backward_translation.size());
                    dej_assert(v_from >= 0);
                    dej_assert((unsigned int)v_to < backward_translation.size());
                    dej_assert(v_to >= 0);
                    dej_assert(orig_v_from < domain_size);
                    dej_assert(orig_v_from >= 0);
                    dej_assert(orig_v_to < domain_size);
                    dej_assert(orig_v_to >= 0);
                    dej_assert(automorphism[orig_v_from] == orig_v_from);
 
                    const int to_recovery_string_start_pt = baked_recovery_string_pt[v_to].first;
                    const int to_recovery_string_end_pt   = baked_recovery_string_pt[v_to].second;
                    const int to_recovery_string_size     = to_recovery_string_end_pt - to_recovery_string_start_pt;
                    const int from_recovery_string_start_pt = baked_recovery_string_pt[v_from].first;
 
                    if(to_recovery_string_size == 0 || baked_recovery_string[to_recovery_string_start_pt] != INT32_MAX) {
                        automorphism[orig_v_from] = orig_v_to;
                        automorphism_supp.push_back(orig_v_from);
                    }
 
                    dej_assert((baked_recovery_string_pt[v_from].second - from_recovery_string_start_pt) ==
                                to_recovery_string_size);
 
                    for (int j = 0; j < to_recovery_string_size; ++j) {
                        const int v_from_t = baked_recovery_string[from_recovery_string_start_pt + j];
                        const int v_to_t   = baked_recovery_string[to_recovery_string_start_pt + j];
 
                        if(v_from_t == INT32_MAX) continue;
 
                        dej_assert(v_from_t >  0?v_to_t >= 0:true);
                        dej_assert(v_to_t   >  0?v_from_t >= 0:true);
                        dej_assert(v_from_t <  0?v_to_t <= 0:true);
                        dej_assert(v_to_t   <  0?v_from_t <= 0:true);
                        if(v_from_t >= 0 && v_to_t >= 0) {
                            dej_assert(automorphism[v_from_t] == v_from_t);
                            automorphism[v_from_t] = v_to_t;
                            automorphism_supp.push_back(v_from_t);
                        } else {
                            const int abs_v_from_t = abs(v_from_t);
                            const int abs_v_to_t   = abs(v_to_t);
 
                            aux_automorphism[abs_v_from_t] = abs_v_to_t;
                            aux_automorphism_supp.push_back(abs_v_from_t);
 
                            use_aux_auto = true;
                        }
                    }
                }
            }
 
            if(use_aux_auto) {
                for(int i = 0; i < aux_automorphism_supp.cur_pos; ++i) {
                    const int v_from = aux_automorphism_supp[i];
                    before_move[v_from]  = automorphism[aux_automorphism[v_from]];
                }
                for(int i = 0; i < aux_automorphism_supp.cur_pos; ++i) {
                    const int v_from = aux_automorphism_supp[i];
                    if(automorphism[v_from] == v_from)
                        automorphism_supp.push_back(v_from);
                    automorphism[v_from] = before_move[v_from];
                }
                reset_automorphism(aux_automorphism.get_array(), aux_automorphism_supp.cur_pos,
                                   aux_automorphism_supp.get_array());
                aux_automorphism_supp.reset();
            }
 
            const int end_pos = automorphism_supp.cur_pos;
 
            if(edge_attached_information > 0) {
                for (int i = 0; i < end_pos; ++i) {
                    recovery_attached_edges_to_automorphism(automorphism_supp[i], automorphism.get_array(),
                                                            &automorphism_supp, &aux_automorphism, &before_move);
                }
            }
 
            if(hook != nullptr) {
                (*hook)(domain_size, automorphism.get_array(), automorphism_supp.cur_pos, automorphism_supp.get_array());
            }
            reset_automorphism(automorphism.get_array(), automorphism_supp.cur_pos, automorphism_supp.get_array());
            automorphism_supp.reset();
        }
 
    private:
 
        // initialize some data structures
        void assure_ir_quotient_init(dejavu::sgraph *g) {
            if (!ir_quotient_component_init) {
                touched_color_cache.initialize(g->v_size);
                ir_quotient_component_init = true;
            }
        }
 
        // deletes edges connected to discrete vertices, and marks discrete vertices for deletion later
        void del_discrete_edges_inplace(dejavu::sgraph *g, coloring *col) {
            int rem_edges = 0;
            int discrete_vert = 0;
            del.reset();
            for (int i = 0; i < col->domain_size;) {
                const int col_sz = col->ptn[i];
                if (col_sz == 0) {
                    ++discrete_vert;
                    del.set(col->lab[i]);
                }
                i += col_sz + 1;
            }
 
            for (int v = 0; v < g->v_size; ++v) {
                if (del.get(v)) {
                    dej_assert(col->ptn[col->vertex_to_col[v]] == 0);
                    rem_edges += g->d[v];
                    g->d[v] = 0;
                    continue;
                }
 
                for (int n = g->v[v]; n < g->v[v] + g->d[v];) {
                    const int neigh = g->e[n];
                    dej_assert(neigh >= 0 && neigh < g->v_size);
                    if (del.get(neigh)) { // neigh == -1
                        const int swap_neigh = g->e[g->v[v] + g->d[v] - 1];
                        //g->e[g->v[v] + g->d[v] - 1] = neigh; // removed this operation because unnecessary?
                        g->e[n] = swap_neigh;
                        --g->d[v];
                        ++rem_edges;
                    } else {
                        ++n;
                    }
                }
            }
            dej_assert(rem_edges % 2 == 0);
        }
 
        void order_according_to_color(dejavu::sgraph *g, int* colmap) {
            bool in_order = true;
            for(int i = 0; i < g->v_size-1 && in_order; ++i) in_order = in_order && (colmap[i] <= colmap[i+1]);
            if(in_order) {
                order_edgelist(g);
                return;
            }
 
            g->initialize_coloring(&c, colmap);
            dejavu::worklist old_arr(g->v_size);
 
            std::memcpy(old_arr.get_array(), g->v, g->v_size*sizeof(int));
            for(int j = 0; j < g->v_size; ++j) {
                g->v[j] = old_arr[c.lab[j]];
            }
 
            std::memcpy(old_arr.get_array(), g->d, g->v_size*sizeof(int));
            for(int j = 0; j < g->v_size; ++j) {
                old_arr[j] = g->d[j];
            }
            for(int j = 0; j < g->v_size; ++j) {
                g->d[j] = old_arr[c.lab[j]];
            }
 
 
            for(int i = 0; i < g->v_size; ++i) {
                colmap[i] = c.vertex_to_col[c.lab[i]];
            }
 
            for(int i = 0; i < g->v_size; ++i) {
                const int map_to = c.lab[i];
                old_arr[map_to] = i; // iso^-1
            }
 
            for(int j = 0; j < g->e_size; ++j) {
                g->e[j] = old_arr[g->e[j]];
            }
 
            dej_assert((int)backward_translation_layers[backward_translation_layers.size() - 1].size() == g->v_size);
            for(int i = 0; i < g->v_size; ++i) {
                old_arr[i] = backward_translation_layers[backward_translation_layers.size() - 1][i];
            }
 
            for(int i = 0; i < g->v_size; ++i) {
                backward_translation_layers[backward_translation_layers.size() - 1][i] = old_arr[c.lab[i]];
            }
 
            order_edgelist(g);
        }
 
    public:
 
        void reduce(dejavu::static_graph *g, dejavu_hook* hook, std::vector<preop> *schedule = nullptr) {
            reduce(g->get_sgraph(), g->get_coloring(), hook, schedule);
        }
 
        void reduce(dejavu::sgraph *g, int *colmap, dejavu_hook* hook, const std::vector<preop> *schedule = nullptr) {
            const std::vector<preop> default_schedule =
                    {deg01, qcedgeflip, deg01, twins, deg2ma, deg2ue, densify2}; // deg2ma, deg2ue, densify2
            // deg01, qcedgeflip, deg2ma, deg2ue, redloop, densify2
            //deg01, qcedgeflip, deg2ma, deg2ue, probe2qc, deg2ma, probeqc, deg2ma, redloop, densify2
 
            if(schedule == nullptr) schedule = &default_schedule;
            if(print) print->print_header();
 
            domain_size = g->v_size;
            saved_hook = hook;
            save_preprocessor = this;
            if(g->v_size == 0)
                return;
            g->dense = !(g->e_size < g->v_size || g->e_size / g->v_size < g->v_size / (g->e_size / g->v_size));
 
            const int test_d = g->d[0];
            int k;
            for (k = 0; k < (g->v_size) && (g->d[k] == test_d); ++k);
 
            int test_col = colmap[0];
            int l;
            for (l = 1; l < (g->v_size) && (colmap[l] == test_col); ++l);
 
            if(k == g->v_size && l == g->v_size) {
                // graph is regular
                skipped_preprocessing = true;
                if(print) print->timer_print("regular", g->v_size, g->e_size);
                return;
            }
 
            backward_translation_layers.emplace_back();
            const size_t back_ind = backward_translation_layers.size() - 1;
            translation_layers.emplace_back();
            backward_translation_layers[back_ind].reserve(g->v_size);
            for (int i = 0; i < g->v_size; ++i) backward_translation_layers[back_ind].push_back(i);
            edge_scratch.allocate(g->e_size);
 
            if(g->v_size > 64)   order_according_to_color(g, colmap);
            else                 order_edgelist(g);
            if(print) print->timer_print("order", g->v_size, g->e_size);
            g->initialize_coloring(&c, colmap);
 
            const int pre_v_size = g->v_size;
            const int pre_e_size = g->e_size;
            const int pre_cells  = c.cells;
 
            dejavu::ir::refinement R_stack;
            if(R1 == nullptr) R1 = &R_stack;
 
            R1->refine_coloring_first(g, &c, -1);
            const bool color_refinement_effective = pre_cells != c.cells;
 
            if (c.cells == g->v_size) {
                //PRINT("(prep-red) graph is discrete");
                g->v_size = 0;
                g->e_size = 0;
                if(print) print->timer_print("discrete", g->v_size, g->e_size);
                return;
            }
 
            add_edge_buff.clear();
            add_edge_buff.reserve(domain_size);
            for (int i = 0; i < domain_size; ++i)
                add_edge_buff.emplace_back(); // do this smarter... i know how many edges end up here
            add_edge_buff_act.initialize(domain_size);
 
            translate_layer_fwd.reserve(g->v_size);
 
            automorphism.allocate(g->v_size);
            for (int i = 0; i < g->v_size; ++i) {
                automorphism.push_back(i);
            }
            automorphism_supp.allocate(g->v_size);
            aux_automorphism.allocate(g->v_size);
            for (int i = 0; i < g->v_size; ++i) {
                aux_automorphism.push_back(i);
            }
            aux_automorphism_supp.allocate(g->v_size);
            _automorphism.allocate(g->v_size);
            _automorphism_supp.allocate(g->v_size);
            for (int i = 0; i < g->v_size; ++i)
                _automorphism[i] = i;
 
            before_move.allocate(domain_size);
 
            worklist_deg0.allocate(g->v_size);
            worklist_deg1.allocate(g->v_size);
 
            domain_size = g->v_size;
 
            // assumes colmap is array of length g->v_size
            del = dejavu::markset();
            del.initialize(g->v_size);
 
            recovery_strings.reserve(g->v_size);
            for (int j = 0; j < domain_size; ++j) recovery_strings.emplace_back();
 
            // eliminate degree 1 + 0 and discrete vertices
            del_discrete_edges_inplace(g, &c);
 
            bool has_deg_0 = false;
            bool has_deg_1 = false;
            bool has_deg_2 = false;
            bool has_discrete = false;
            bool graph_changed = false;
            int count_deg0 = 0;
            int count_deg1 = 0;
            int count_deg2 = 0;
            int count_discrete = 0;
 
            for (int i = 0; i < c.domain_size;) {
                const int v = c.lab[i];
                const int col_sz = c.ptn[i] + 1;
                switch (g->d[v]) {
                    case 0:
                        count_deg0 += col_sz;
                        break;
                    case 1:
                        count_deg1 += col_sz;
                        break;
                    case 2:
                        count_deg2 += col_sz;
                        break;
                    default:
                        break;
                }
                count_discrete += (col_sz == 1);
                i += col_sz;
            }
            copy_coloring_to_colmap(&c, colmap);
            if(print) print->timer_print("colorref", g->v_size, c.cells);
 
            //if(!has_deg_0 && !has_deg_1 && !has_deg_2 && !has_discrete) return;
            has_deg_0 = (count_deg0 > 4);
            has_deg_1 = (count_deg1 > 8);
            has_deg_2 = (count_deg2 > 128); /*< if there's only very few, there's really no point... */
            has_discrete = (count_discrete > 4);
 
            if(!has_deg_0 && !has_deg_1 && !has_deg_2 && !has_discrete) return;
 
            if (schedule != nullptr) {
                del_e = dejavu::markset();
                del_e.initialize(g->e_size);
 
                for (size_t pc = 0; pc < schedule->size(); ++pc) {
                    if (g->v_size <= 1) {
                        return;
                    }
                    preop next_op = (*schedule)[pc];
                    const int pre_v = g->v_size;
                    const int pre_e = g->e_size;
                    switch (next_op) {
                        case preop::deg01: {
                            if(!has_deg_0 && !has_deg_1 && !has_discrete && !graph_changed) break;
                            red_deg10_assume_cref(g, colmap, hook);
                            perform_del(g, colmap);
                            if(print) print->timer_print("deg01", g->v_size, g->e_size);
 
                            if(!(g->v_size == pre_v_size && g->e_size == pre_e_size && !color_refinement_effective)) {
                                order_according_to_color(g, colmap);
                                if(print) print->timer_print("order", g->v_size, g->e_size);
                            }
                            dej_assert(_automorphism_supp.cur_pos == 0);
                            break;
                        }
                        case preop::deg2ma: {
                            if(!has_deg_2 && !graph_changed) break;
                            red_deg2_path_size_1(g, colmap);
                            perform_del_add_edge(g, colmap);
                            if(print) print->timer_print("deg2ma", g->v_size, g->e_size);
                            dej_assert(_automorphism_supp.cur_pos == 0);
 
                            break;
                        }
 
                        case preop::twins: {
                            // check for twins and mark them for deletion
                            //int s_twins = remove_twins(g, colmap, hook);
                            const int s_twins_true  = twins_partition_refinement(g, colmap, hook, false);
                            if(s_twins_true > 0) perform_del(g, colmap);
                            const int s_twins_false = twins_partition_refinement(g, colmap, hook, true);
                            if(s_twins_false > 0) perform_del(g, colmap);
                            const int s_twins = s_twins_true + s_twins_false;
 
                            if(s_twins > 0) {
                                // might distinguish vertices, so let's apply color refinement again
                                g->initialize_coloring(&c, colmap);
                                R_stack.refine_coloring_first(g, &c);
                                copy_coloring_to_colmap(&c, colmap);
 
                                dej_assert(_automorphism_supp.cur_pos == 0);
 
                                // twins can re-activate other routines, so let's start over...
                                if(s_twins >= 16)  pc = 0;
                            }
 
                            if(print) print->timer_print("twins", g->v_size, g->e_size);
                            break;
                        }
                        case preop::deg2ue: {
                            if(!has_deg_2 && !graph_changed) break;
 
                            has_discrete = red_deg2_color_cycles(g, colmap);
 
                            red_deg2_unique_endpoint(g, colmap);
                            perform_del_add_edge(g, colmap);
 
                            red_deg2_trivial_connect(g, colmap);
                            perform_del_add_edge(g, colmap);
                            if(has_discrete) {
                                perform_del_discrete(g, colmap);
                            }
 
                            if(print) print->timer_print("deg2ue", g->v_size, g->e_size);
                            dej_assert(_automorphism_supp.cur_pos == 0);
                            break;
                        }
                        case preop::densify2: {
                            if(!has_deg_2 && !graph_changed) break;
                            red_deg2_densify(g, colmap);
                            perform_del_add_edge_general(g, colmap);
 
                            red_deg2_densifysc1(g, colmap);
                            perform_del_add_edge_general(g, colmap);
 
                            if(print) print->timer_print("densify2", g->v_size, g->e_size);
                            dej_assert(_automorphism_supp.cur_pos == 0);
                            break;
                        }
                        case preop::qcedgeflip: {
                            red_quotient_edge_flip(g, colmap);
                            perform_del_edge(g);
 
                            if(print) print->timer_print("qcedgeflip", g->v_size, g->e_size);
                            break;
                        }
                    }
                    graph_changed = graph_changed || pre_v != g->v_size || pre_e != g->e_size;
                }
            }
 
            skipped_preprocessing = g->v_size == domain_size;
        }
 
        void save_my_hook(dejavu_hook *hook) {
            saved_hook = hook;
        }
 
        // bliss usage specific:
#if defined(BLISS_VERSION_MAJOR) && defined(BLISS_VERSION_MINOR)
#if ( BLISS_VERSION_MAJOR >= 1 || BLISS_VERSION_MINOR >= 76 )
        void bliss_hook(unsigned int n, const unsigned int *aut) {
          auto p = save_preprocessor;
          p->pre_hook_buffered(n, (const int *) aut, -1, nullptr, p->saved_hook);
       }
#else
        static inline void bliss_hook(void *user_param, unsigned int n, const unsigned int *aut) {
                    auto p = (preprocessor *) user_param;
                    p->pre_hook_buffered(n, (const int *) aut, -1, nullptr, p->saved_hook);
                }
#endif
#else
        static inline void bliss_hook(void *user_param, unsigned int n, const unsigned int *aut) {
                auto p = (preprocessor *) user_param;
                p->pre_hook_buffered(n, (const int *) aut, -1, nullptr, p->saved_hook);
            }
#endif
        // Traces usage specific:
        static inline void traces_hook(int, int* aut, int n) {
            auto p = save_preprocessor;
            p->pre_hook_buffered(n, (const int *) aut, -1, nullptr, p->saved_hook);
        }
 
        void traces_save_my_preprocessor() {
            save_preprocessor = this;
        }
 
        // nauty usage specific:
        static inline void nauty_hook(int, int* aut, int*, int, int, int n) {
            auto p = save_preprocessor;
            p->pre_hook_buffered(n, (const int *) aut, -1, nullptr, p->saved_hook);
        }
 
        void nauty_save_my_preprocessor() {
            save_preprocessor = this;
        }
 
        // saucy usage specific:
        static inline int saucy_hook(int n, const int* aut, int nsupp, int* supp, void* user_param) {
            auto p = (preprocessor *) user_param;
            p->pre_hook_buffered(n, (const int *) aut, nsupp, supp, p->saved_hook);
            return true;
        }
 
        // dejavu usage specific
        static inline void _dejavu_hook(int n, const int* aut, int nsupp, const int* supp) {
            auto p = save_preprocessor;
            if(p->skipped_preprocessing && !p->decomposer) {
                if(p->saved_hook != nullptr) {
                    (*p->saved_hook)(n, aut, nsupp, supp);
                }
                return;
            }
            p->pre_hook_buffered(n, (const int *) aut, nsupp, supp, p->saved_hook);
        }
    };
}
#endif //SASSY_H