decompose.c

//
// Created by Stephane on 10/03/2020.
//

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "decompose.h"
#include "components.h"
#include "compression.h"
#include "swaps.h"
#include "graph.h"
#include "main.h"
#include "separator.h"
#include "tree.h"
#include "utils.h"
#include "heap.h"

extern int trace;
extern char * last_name;


int *bestTree;
int bestHeight;
int bestEvalMode = NONE;


Separator separator;

Node * theNodes;

SET * theSets;
int nbAllocatedSets;

int ** components;
int ** componentsOf;
int nbAllocatedComponents;
int topComponents;



int nbCallsDecompose = 0;
int nbCallsDecomposeComponents = 0;

extern int indpdNbCalls;

int mustVerifyDecomposition = 0;



void updateBestDecomposition(Node root, Graph g) {
    if ((root != NULL) && (root->height < bestHeight)) {
        bestHeight = root->height;
        saveToTable(bestTree, g);
        bestEvalMode = modeEvalSeparator;
    }
}



void allocAndInitializeDecomposition(Graph g) {
    separator = newSeparator(g->n, g);

    allocSeparation(g);
    allocNodes(g);
    nbAllocatedSets = 100;
    allocSets(g, nbAllocatedSets);

    nbAllocatedComponents = 100;
    allocComponents(g, nbAllocatedComponents);
    allocSearchConnectedComponents(g);

    bestTree = malloc(g->n * sizeof(int));
    bestHeight = g->n+1;
}

void initializeRun(Graph g, int S[], int pos[]) {
    topComponents = 0;

    for (int i = 0; i < g->n; i ++) {
        S[i] = i;
        pos[i] = i;
    }
}

void markRemoved(int S[], int n, int pos[]) {
    for (int i = 0; i < n; i ++)
        pos[S[i]] = NONE;
}


#ifdef STATS_SEPARATION
extern clock_t *timesSep;
extern int *nbSearches;
extern int *sumSizes;
#endif


void testDecompose(Graph g, int nbRuns, int algo) {
    int *S = malloc(g->n * sizeof(int)); // list of vertices
    int *posS = malloc(g->n * sizeof(int)); // booleans for separation
    int *nbCriticalNodes = malloc(g->n * sizeof(int)); // to calculate the numbers of critical nodes at each level
    int bestHDec = g->n+1;
    Node root;
    int h0, h1, h2;
    int nbEffectiveRuns = 0;
    int nbRunsAttempted = 0;
    int trace = 0;

    if (perForceChoiceNotInC == NONE)
        perForceChoiceNotInC = (g->n > 7600) ? 10 : 0;

    clock_t startTimeRun;
    clock_t startTimeImprove;
    clock_t totalTimeDecompose = 0;
    clock_t totalTimeImprove = 0;

    allocAndInitializeDecomposition(g);
    allocConnectionsHeap(g);

    if (preliminaryGreedyDecomposition && (algo != GREEDY_ALGORITHM) && (g->n < maxSizePrelimGreedy)) {
        nbEffectiveRuns ++;
        initConHeap(g);
        root = greedyDecompose(conHeap, 0, g, g->n+1, 1);
        if (trace) printf("greedy first decomp h=%d  %.1fs\n", root->height, (float) (clock()-startTime)/CLOCKS_PER_SEC);
        if (root != NULL) updateBestDecomposition(root, g);
        if (stopSearch) goto SAVE_AND_EXIT;
    }

    for (int i = 0; i < nbRuns; i ++) {

        startTimeRun = clock();
        if ((startTimeRun-startTime)/CLOCKS_PER_SEC >= time_limit)
            break;

        if (algo == SEPARE_AND_EXPLORE) {
            initializePriorities(g);
            initializeRun(g, S, posS);
        }

        int maxDepth = bestHDec*(100+depthThreshold)/100;
        if (noSwapsInFirstRuns) maxDepth = bestHDec+1;

        initConHeap(g);

        if (algo == GREEDY_ALGORITHM) {
            root = greedyDecompose(conHeap, 0, g, maxDepth, 1);
            if (trace) printf("greedy decomp h=%d  %.1fs\n", (root == NULL) ? -1 : root->height, (float) (clock()-startTime)/CLOCKS_PER_SEC);
        }
        else if (algo == SEPARE_AND_EXPLORE) {
            root = decompose(S, g->n, posS, 0, g, maxDepth, 0);
            if (trace) printf("decompose h=%d  %.1fs\n", (root == NULL) ? -1 : root->height, (float) (clock()-startTime)/CLOCKS_PER_SEC);
        }

        nbRunsAttempted ++;
        if (stopSearch) break;

        totalTimeDecompose += (clock()-startTimeRun);

        if (root == NULL) {
            if (nbRunsSeparation < 50) nbRunsSeparation ++;
            continue;
        }

        nbEffectiveRuns ++;
        h0 = root->height;
        if (root->height < bestHDec) bestHDec = root->height;

        if (printStats) printStatsOnCriticalNodes(root, nbCriticalNodes);

        updateBestDecomposition(root, g);

        if (noSwapsInFirstRuns) {
            if ((clock()-startTime) / CLOCKS_PER_SEC < timeWithNoSwaps * time_limit / 100) continue; // Shunt improve at the beginning
            root = buildTreeFromTable(bestTree, g);
            noSwapsInFirstRuns = 0;
            nbRuns = 1000000;
        }

        startTimeImprove = clock();
        if (pullUpIndependentSubtrees || searchSwapsInCriticalBranch) {
            root = makeSwapsInCriticalBranch(root, g);
            if (trace) printf("improve  %d --> %d  %.1fs\n", h0, root->height, (float) (clock()-startTime)/CLOCKS_PER_SEC);
            updateBestDecomposition(root, g); // useless
        }
        totalTimeImprove += (clock()-startTimeImprove);

        if (stopSearch) break;

        if (mustVerifyDecomposition) if ( ( ! verifyTree(root, 0)) || ( ! verifyDecomposition(root, g)) ) exit(0);

        if (0) printf("decompose h=%d --> %d  %.1fs \n", h0, root->height, (float) (clock()-startTimeRun)/CLOCKS_PER_SEC);

#ifdef THE_OLD_IMPROVEMENT
        if ( ! improve) { if (trace) printf("\n"); continue; }
        if (root->height > bestHDec*(100+depthThreshold)/100) {if (trace) printf("abandon\n"); continue;}
    IMPROVE_TREE:

        if (stopSearch) break;

        if (compressTree) {
            exploreAndCompress(root, g, 0);
            h1 = root->height;
            verifyTree(root, 0);

            if (trace) printf("h1=%d %.1fs %d  ", h1, (float) (clock() - startTimeRun) / CLOCKS_PER_SEC, indpdNbCalls);
            if (root->height < bestHeight) {
                bestHeight = root->height;
                saveToTable(bestTree, g);
                if (trace) printf("\nbestC=%d ", root->height);
            }

            if (stopSearch) break;
        }

        searchForksSwaps(root, g, root->height);
        h2 = root->height;
        verifyTree( root, 0);

        if (trace) printf("h2=%d %.1fs %d  ", h2, (float) (clock()-startTimeRun)/CLOCKS_PER_SEC, indpdNbCalls);
        if (root->height < bestHeight) { bestHeight = root->height; saveToTable(bestTree, g); if (trace) printf("\nbestFS=%d ", root->height); }

        if (h2 < h0) goto IMPROVE_TREE;

        if (trace) printf("\n");
#endif
    }

#ifdef STATS_SEPARATION
    for (int i = 0; ; i ++) {
        if (nbSearches[i] == 0) break;
        printf("[%d]  %d x %.1fs   (size=%d)\n", i, nbSearches[i],
                (float) (timesSep[i] / nbSearches[i]) / CLOCKS_PER_SEC, sumSizes[i]/nbSearches[i]);
    }
#endif

SAVE_AND_EXIT:
    if (printResult) printf("%s %d %d   runs=%d/%d   calls=%d+%d+%d   h=%d    %.2fs + %.2fs   hdec=%d  mode=%d\n",
            last_name, g->n, g->m, nbEffectiveRuns, nbRunsAttempted,
            nbCallsDecompose, nbCallsDecomposeComponents, nbCallsSwapNode,
            bestHeight,
            (float) totalTimeDecompose/CLOCKS_PER_SEC, (float) totalTimeImprove/CLOCKS_PER_SEC,
            bestHDec, bestEvalMode);
    if (printSolution) {
        PACEOutput(fileSolution, bestHeight, bestTree, g);
    }
}


//
// Decomposition : SEPARE_AND_EXPLORE
//

Node decompose(int *S, int n, int pos[], int depth, Graph g, int hmax, int connected) {

    nbCallsDecompose ++;
    //int nbCalls = nbCallsDecompose;

    if (trace) printf("[%d] call %d  n=%d   \n", depth, nbCallsDecompose, n);

    assert(depth <= g->n);

    if (stopSearch) return NULL;
    if (n == 0) return NULL;
    if (hmax <= 0) return NULL;

    if (n <= 3) { markRemoved(S, n, pos); return makeSmallTree(S, n, g); }

    if (0 && n <= sizeSwitchToGreedy) {
        conHeap->val = S;
        conHeap->n = n;
        // must first rebuild conHeap
        return greedyDecompose(conHeap, depth, g, hmax, 1);
    }

    SET set = NULL;
    if (n > 105) {
        if (theSets[depth] == NULL) reAllocSets(100, g);
        set = makeSet(S, n, theSets[depth]);
    }
    else
        Introsort(S, S, S + n - 1);

    // Search connected components
    if (depth > 0) {
        int nbComp = searchConnectedComponents(set, S, n, g);
        if (stopSearch) return NULL;
        if (nbComp > 1) {
            if (trace) printf("components: %d   ", nbComp);
            return decomposeConnectedComponents(iComp, nbComp, S, n, pos, depth, g, hmax);
        }
    }

    int nbEdges = initializeNbNeighbors(set, S, n, pos, g); // for searchSeparator()

    if (2*nbEdges == n*(n-1)) {
        markRemoved(S, n, pos);
        return makeSingleBranch(S, n, NULL, g);
    }

    //searchSeparator(S, n, set, g, separator, (n < 20) ? 1 : nbRunsSeparation, (n < 20) ? 1 : nbFlushes, NULL, 0);
    searchSeparator(S, n, set, g, separator, nbRunsSeparation, nbFlushes, depth);

    if (stopSearch) return NULL;

    if (0) {
        improveSeparation(separator, g, 0); // Seg fault with 100 runs "./treedepth -file ../public/heur_007.gr "
        printf("improve :: %d %d %d\n", separator->B->n, separator->C->n, separator->A->n);
    }

    //root = makeSingleBranchWithDisconnectedVertices(separator->C->val, separator->C->n, separator->nbABDV, g);
    Node root = makeSingleBranch(separator->C->val, separator->C->n, nbNInABCopy, g);

    Node theLastOfTheBranch = bottomNode;
    markRemoved(separator->C->val, separator->C->n, pos);

    int nA = separator->A->n;
    int nB = separator->B->n;
    int nC = separator->C->n;
    Node A, B;

    if (0) { for (int i = 0; i < depth; i ++) printf(" "); printf("%d %d     %d\n", 1000*nA/n, 1000*nB/n, 1000*nC/n); }

    assert(nC > 0); // Indeed, the subgraph is connected here.


    int *listA = S;
    int *listB = S+nA; // no need for S after here

    if (trace) printf("   sep: %d,%d,%d\n", nA, nB, nC);

    copyList(listA, separator->A->val, nA);
    copyList(listB, separator->B->val, nB);

    // Here nbNeighborsInA and nbNeighborsInB contain for each vertex in listA and listB
    // its number of neighbors in listA and listB. No need to recalculate in initSeparator.

    if (nA >= nB) {
        A = decompose(listA, nA, pos, depth + 1, g, hmax - nC, 0);
        if (stopSearch) return NULL;
        if ((A == NULL) && (nA > 0)) {
            root = NULL;
            goto FIN;
        }
        B = decompose(listB, nB, pos, depth + 1, g, hmax - nC, 0);
        if (stopSearch) return NULL;
        if ((B == NULL) && (nB > 0)) {
            root = NULL;
            goto FIN;
        }
    }
    else {
        B = decompose(listB, nB, pos, depth + 1, g, hmax - nC, 0);
        if (stopSearch) return NULL;
        if ((B == NULL) && (nB > 0)) {
            root = NULL;
            goto FIN;
        }
        A = decompose(listA, nA, pos, depth + 1, g, hmax - nC, 0);
        if (stopSearch) return NULL;
        if ((A == NULL) && (nA > 0)) {
            root = NULL;
            goto FIN;
        }
    }


    if (0) if (sizeForest(A)+ sizeForest(B)+ sizeForest(root) != nA+nB+nC) {
        printf("nA=%d nB=%d nC=%d  A:%d  B:%d root:%d\n", nA, nB, nC, sizeForest(A), sizeForest(B), sizeForest(root));
        exit(0);
    }


    if (root == NULL) {
        if (A == NULL) { root = B; goto FIN; }
        if (B == NULL) { root = A; goto FIN; }
        //assert(A->next == NULL);
        //assert(B->next == NULL);
        root = A;
        while (A->next != NULL) A = A->next;
        A->next = B;
        assert((root->next == NULL) || (depth > 0));
        goto FIN;
    }


    if (A != NULL) {
        if (A->next == NULL)
            addChild(A, theLastOfTheBranch);
        else {
            while (A != NULL) {
                Node next = A->next;
                addChild(A, theLastOfTheBranch);
                A = next;
            }
        }
    }


    if (B != NULL) {
        if (B->next == NULL)
            addChild(B, theLastOfTheBranch);
        else {
            while (B != NULL) {
                Node next = B->next;
                addChild(B, theLastOfTheBranch);
                B = next;
            }
        }
    }

    if (theLastOfTheBranch->father != NULL)
        updateIncHeightOnBranch(theLastOfTheBranch->father, theLastOfTheBranch->height);

    //updateBranchUp(theLastOfTheBranch, root);
    FIN:

    //verifyTree(root, depth);
    return root;
}


Node decomposeConnectedComponents(int *iComp, int nbComp, int *S, int n, int pos[], int depth, Graph g, int hmax) {

    Node list = NULL, node;

    nbCallsDecomposeComponents ++;

    int *sizes = malloc(nbComp*sizeof(int)); // To memorize sizes after first explorations

    sortListByComponent(S, n, nbComp, sizes);
    int nbV = 0;

    for (int i = 0; i < nbComp; i ++) {
        if (stopSearch) { list = NULL; break; }

        node = decompose(S+nbV, sizes[i], pos, depth, g, hmax, 1);

        if ((node == NULL) || (stopSearch)) return NULL;

        node->next = list;
        if (list != NULL) list->prev = node;
        list = node;

        nbV += sizes[i];
    }

    return list;

#ifdef OLD_VERSION
    int *component, *componentOf;

    if (components[topComponents] == NULL) reAllocComponents(g, 100);
    component = components[topComponents];
    componentOf = componentsOf[topComponents];
    topComponents ++;

    // Memorize iComp[] that may be used in next calls
    for (int i = 0; i < n; i ++)
        componentOf[i] = iComp[S[i]];

    for (int i = 0; i < nbComp; i ++) {
        if (stopSearch) { list = NULL; break; }

        int size = extractList(i, S, n, componentOf, component);

        node = decompose(component, size, depth, g, hmax, 1);
        if (node == NULL) { topComponents --; return NULL; }

        node->next = list;
        if (list != NULL) list->prev = node;
        list = node;
    }

    topComponents --;
    return list;
#endif
}



//
// Greedy-decompose (V,E)
// 1. select vertex u with highest degree,
// 2. C:=disconnectedComponents(V-{u},E),
// 3. tree:=<u>,
// 4. for each c in C do
// 5.    addChild(Greedy-decompose(c,E), tree),
// 6. return tree,

int lastRemovedVertex = NONE;

Node greedyDecompose(Heap conHeap, int depth, Graph g, int hmax, int connected) {

    nbCallsDecompose ++;

    if (depth > maxDepthGreedy) return NULL;

    //int localNbCalls = nbCallsDecompose;

    if (trace) printf("[%d] call %d  n=%d  max=%d ", depth, nbCallsDecompose, conHeap->n, maxNbN4Con);

    if (stopSearch) return NULL;
    if (conHeap->n == 0) return NULL;
    if (hmax == 0) return NULL;

    if (conHeap->n <= SIZE_SMALL_TREE) {
        Node root = makeSmallTree(conHeap->val, conHeap->n, g);
        heapRemoveAll(conHeap);
        return root;
    }

    // Isolated vertices
    if (nbUnconnectedVertices > 0) {
        //assert(connected == 0);
        if (trace) printf(" :: %d unconnected vertices :: %d --> %d\n", nbUnconnectedVertices, conHeap->n, conHeap->n - nbUnconnectedVertices);
        Node last;
        Node root =  makeListWithUnconnectedVertices(conHeap, g, &last);
        if (conHeap->n > 0) {
            nbUnconnectedVertices = 0;
            Node sibling = greedyDecompose(conHeap, depth + 1, g, hmax, 0);
            if (sibling == NULL)
                return NULL;
            last->next = sibling;
        }
        return root;
    }

    // Search disconnected components
    if (( ! connected) && (depth > 0) && (lastRemovedVertex != NONE) &&
            (nbN4Con[lastRemovedVertex] != 1)) { // added 7/04/2020
        int nbComp = searchConnectedComponentsGreedy(NULL, conHeap, g, lastRemovedVertex);
        if (nbComp > 1) {       //if (trace) { printf("components [%d]   ", nbComp); for (int i = 0; i < nbComp; i ++) printf("%d,", compSizes[i]); printf("\n");}
            return greedyDecomposeComponents(iComp, nbComp, conHeap, depth, g, hmax);
        }
    }

    if (stopSearch) return NULL;

    // Build a tree whose root is a node with highest degree

    int vertex = vertexWithManyConnections(conHeap, MAX_CONNECTED_BEST_AT_RANDOM);

    if (trace) printf("remove %d (nbN=%d)\n", vertex, nbN4Con[vertex]);

    //if (nbInListWithThisValue(maxNbN4Con, conHeap->val, conHeap->n, nbN4Con) != nbMax4Con) { printf("AAA max=%d  %d et non pas %d\n",maxNbN4Con, nbInListWithThisValue(maxNbN4Con, conHeap->val, conHeap->n, nbN4Con), nbMax4Con); exit(0); }

    heapRemove(vertex, conHeap);

    if (nbN4Con[vertex] == maxNbN4Con) {
        if (--nbMax4Con == 0) updateNbMaxConInHeap(conHeap);
    }

    if (nbN4Con[vertex] == 0) nbUnconnectedVertices --;

    //printf("1. max=%d nbMax=%d  -->  ", maxNbN4Con, nbMax4Con);
    //if (nbInListWithThisValue(maxNbN4Con, conHeap->val, conHeap->n, nbN4Con) != nbMax4Con) { printf("BBB max=%d  %d et non pas %d\n",maxNbN4Con, nbInListWithThisValue(maxNbN4Con, conHeap->val, conHeap->n, nbN4Con), nbMax4Con); exit(0); }

    updateHeapNbNAfterRemoval(vertex, conHeap, g);

    //assert(conHeap->n > 0);
    if (nbMax4Con == 0) updateNbMaxConInHeap(conHeap);

    //printf("2. max=%d nbMax=%d\n", maxNbN4Con, nbMax4Con);
    //if (nbInListWithThisValue(maxNbN4Con, conHeap->val, conHeap->n, nbN4Con) != nbMax4Con) { printf("CCC max=%d  %d et non pas %d\n",maxNbN4Con, nbInListWithThisValue(maxNbN4Con, conHeap->val, conHeap->n, nbN4Con), nbMax4Con); exit(0); }

    Node root = theNodes[vertex];
    resetNode(root);

    lastRemovedVertex = vertex;

    Node children = greedyDecompose(conHeap, depth+1, g, hmax, 0);

    if (children == NULL)
        return NULL;

    if (children->next == NULL)
        addChild(children, root);
    else {
        while (children != NULL) {
            Node child = children;
            children = children->next;
            addChild(child, root);
        }
    }

    //verifyTree(root, 0);

    return root;
}



Node greedyDecomposeComponents(int iComp[], int nbComp, Heap conHeap, int depth, Graph g, int hmax) {
// Normally there is no component of size 1
    nbCallsDecomposeComponents ++;
    int numCall = nbCallsDecomposeComponents;

    int *sizes = malloc(nbComp*sizeof(int)); // To memorize sizes after first explorations

    sortListByComponent(conHeap->val, conHeap->n, nbComp, sizes);

    if (trace) {
        printf(" :: split call %d : %d components : ", nbCallsDecomposeComponents, nbComp);
        for (int i = 0; i < nbComp; i ++) printf("%d ", sizes[i]);
        printf("\n");
    }

    // The number of neighbors of u in S[] is the number of neighbors of u in its component
    // No need to update.
    Node list = NULL, node;
    int nb = 0;
    int *val = conHeap->val;
    for (int i = 0; i < nbComp; i ++) {

        if (stopSearch) { conHeap->val = val; free(sizes); return NULL; }

        conHeap->n = sizes[i];
        conHeap->val = val+nb;
        nb += sizes[i];

        if (i == nbComp-1) free(sizes);

        lastRemovedVertex = NONE;

        if (conHeap->n <= SIZE_SMALL_TREE) {
            node = makeSmallTree(conHeap->val, conHeap->n, g);
            heapRemoveAll(conHeap);
        }
        else if (rebuildConHeap(conHeap)) {
            // clique
            node = makeSingleBranch(conHeap->val, conHeap->n, NULL, g);
            heapRemoveAll(conHeap);
        }
        else
            node = greedyDecompose(conHeap, depth, g, hmax, 1); // the component is connected

        if (0 && ! verifyDecomposition(node, g)) { printf("Error: decCompCall %d  depth=%d\n", numCall, depth); exit(0); }

        if (node == NULL) { conHeap->val = val; if (i != nbComp-1) free(sizes); return NULL; }

        node->next = list;
        if (list != NULL) list->prev = node;
        list = node;
    }

    conHeap->val = val;
    //free(sizes);

    return list;
}




Node makeListWithUnconnectedVertices(Heap heap, Graph g, Node *last) {
    // unconnected vertices are at the end of the heap
    // They are just removed from the heap, the degree of their already removed neighbors is not updated
    int firstUV = heap->n - nbUnconnectedVertices;
    Node node = NULL;
    *last = NULL;

    for (int i = heap->n - 1; i >= firstUV ; i --) {
        int u = heap->val[i];
        assert (nbN4Con[u] == 0);
        //heapRemove(heap->val[i], heap);
        heap->ind[u] = NONE;
        Node p = theNodes[u];
        resetNode(p);
        if (*last == NULL) *last = p;
        p->next = node;
        node = p;
    }
    heap->n = firstUV;
    return node;
}


Node makeSmallTree(int *S, int n, Graph g) {

    if (n == 0) return NULL;

    Node u = theNodes[S[0]];
    resetNode(u);

    if (n == 1) return u;

    Node v = theNodes[S[1]];
    resetNode(v);

    if (n == 2) {
        if (areNeighbours(u->vertex, v->vertex, g)) {
            addChild(v, u);
            return u;
        }
        //connectSiblings(u, v);
        u->next = v;
        return u;
    }

    Node w = theNodes[S[2]];
    resetNode(w);

    int status = 0;
    if (areNeighbours(u->vertex, v->vertex, g))
        status += 1;
    if (areNeighbours(u->vertex, w->vertex, g))
        status += 2;
    if (areNeighbours(v->vertex, w->vertex, g))
        status += 4;


    if (status == 7) { // Clique -> branch
        addChild(w, v);
        addChild(v, u);
        return u;
    }

    if (status == 6) { // (u,w) (v,w) -> w(u,v)
        addChild(u, w);
        addChild(v, w);
        return w;
    }

    if (status == 4) { // (v,w) -> w(v),u
        addChild(v, w);
        w->next = u; // addChild(u, w);
        return w;
    }

    if (status == 2) { // (u,w) -> w(u),v
        addChild(u, w);
        w->next = v; // addChild(v, w)
        //connectSiblings(w, v);
        return w;
    }

    if (status == 5) { // (u,v) (v,w) -> v(u,w)
        addChild(w, v);
        addChild(u, v);
        return v;
    }

    if (status == 1) { // (u,v) -> v(u),w
        addChild(u, v);
        v->next = w; // addChild(w, v);
        //connectSiblings(v, w);
        return v;
    }

    if (status == 3) { // (u,v) (u,w) -> u(v,w)
        addChild(v, u);
        addChild(w, u);
        return u;
    }

    //addChild(v, u);
    //addChild(w, u);
    //return u;

    // (status == 0) :: (u,v) (u,w) -> u(v,w)
    u->next = v;
    v->next = w;
    //connect3Siblings(u, v, w);

    return u;
}





//
// Utils
//



void allocNodes(Graph g) {
    theNodes = malloc(g->n*sizeof(Node));
    for (int i = 0; i < g->n; i ++)
        theNodes[i] = newNode(i, g);
}

void allocSets(Graph g, int max) {
    theSets = calloc((size_t) (g->n), sizeof(SET));
    for (int i = 0; i < max; i ++)
        theSets[i] = allocSet(g->n + 100);
}


void allocComponents(Graph g, int max) {
    components = calloc((size_t) (g->n), sizeof(int *));
    componentsOf = calloc((size_t) (g->n), sizeof(int *));
    for (int i = 0; i < max; i ++) {
        components[i] = malloc((size_t) (g->n) * sizeof(int));
        componentsOf[i] = malloc((size_t) (g->n) * sizeof(int));
    }
}



// Alloc new sets if necessary
void reAllocSets(int size, Graph g) {

    for (int i = nbAllocatedSets; i < nbAllocatedSets+size; i++) {
        theSets[i] = allocSet(g->n);
    }
    nbAllocatedSets += size;
}


void reAllocComponents(Graph g, int size) {
    for (int i = nbAllocatedComponents; i < nbAllocatedComponents+size; i ++) {
        components[i] = malloc((size_t) (g->n) * sizeof(int));
        componentsOf[i] = malloc((size_t) (g->n) * sizeof(int));
    }
    nbAllocatedComponents += size;
}




void PACEOutput(FILE *file, int height, int fathers[], Graph g) {

    if (file == NULL) file = stdout;

    fprintf(file, "%d\n", height);

    for (int i = 0; i < g->n; i ++) {
        fprintf(file, "%d\n", fathers[i]);
    }
}


void saveToTable(int T[], Graph g) {
    for (int i = 0; i < g->n; i ++) {
        T[i] = (theNodes[i]->father == NULL) ? 0 : (theNodes[i]->father->vertex+1);
    }
}

Node buildTreeFromTable(int T[], Graph g) {
    Node root = NULL;

    for (int i = 0; i < g->n; i ++)
        theNodes[i]->father = theNodes[i]->next = theNodes[i]->fbs = NULL;

    for (int i = 0; i < g->n; i ++) {
        Node father = (T[i] == 0) ? NULL : theNodes[T[i]-1];
        theNodes[i]->father = father;
        if (father != NULL) {
            theNodes[i]->next = father->fbs;
            father->fbs = theNodes[i];
        } else
            root = theNodes[i];
    }
    return root;
}


// Verify that the decomposition tree is correct for the edges of g
// Note : tree may be a partial decomposition of g (does not contain all the vertices of g)
int exploreAndSetFathers(Node p, int F[]) {
    Node q = p->fbs;
    if (F[p->vertex] != NONE) { printf("vertex %d occurs twice !\n", p->vertex); return 0; }
    while (q != NULL) {
        if ( ! exploreAndSetFathers(q, F)) return 0;
        F[q->vertex] = p->vertex;
        q = q->next;
    }
    return 1;
}

int areInSameBranch(int u, int v, int F[]) {
    int i = u;
    while ((i != NONE) && (i != v))
        i = F[i];
    if (i == v) return 1;
    i = v;
    while ((i != NONE) && (i != u))
        i = F[i];
    return (i == u);
}

int verifyDecomposition(Node tree, Graph g) {
    //printf("43e neighbor of %d = %d\n", 13280-1, g->lists[13280-1][43]);
    if (sizeTree(tree) != g->n) { printf("Error size %d nodes\n", sizeTree(tree)); return 0; }
    int F[g->n];
    for (int i = 0; i < g->n; i ++)
        F[i] = NONE;
    if ( ! exploreAndSetFathers(tree, F)) return 0;
    int nbNONE = 0;
    for (int u = 0; u < g->n; u ++) if (F[u] == NONE) nbNONE ++;
    if (nbNONE != 1) { printf("Error : %d nodes not in tree\n", nbNONE-1); return 0; }
        //printf(" in same branch 13280 16624 : %d\n", areInSameBranch(13280-1, 16624-1, F));
    for (int u = 0; u < g->n; u ++) {
        if (F[u] != NONE) {
            for (int *p = g->lists[u]; *p != NONE; p ++) {
                //if ((u == 13280-1) && (*p == 16624-1)) printf("test edge %d %d %d\n", u, *p, (int) (p-g->lists[u]));
                //if ((u == 16624-1) && (*p == 13280-1)) printf("test edge %d %d %d\n", u, *p, (int) (p-g->lists[u]));
                if (F[*p] != NONE) {
                    // g contains edge (u,q->val) and u and q->val are in the tree
                    //if ((u == 16624-1) && (*p == 13280-1)) printf("edge %d %d %d\n", u, *p, (int) (p-g->lists[u]));
                    //if ((u == 13280-1) && (*p == 16624-1)) printf("edge %d %d %d\n", u, *p, (int) (p-g->lists[u]));
                    if ( ! areInSameBranch(u, *p, F)) {
                        printf("Error edge (%d,%d)\n", u, *p);
                        return 0;
                    }
                }
            }
        }
#ifdef NOTDEF
        if (F[u] != NONE) {
            LIST q = g->adj[u];
            while (q != NULL) {
                if ((u == 13280) && (q->val == 16624)) printf("Yahou\n");
                if ((u == 16624) && (q->val == 13280)) printf("Yahou\n");
                if (F[q->val] != NONE) {
                    // g contains edge (u,q->val) and u and q->val are in the tree
                    if ( ! areInSameBranch(u, q->val, F)) {
                        printf("Error edge (%d,%d)\n", u, q->val);
                        return 0;
                    }
                }
                q = q->suiv;
            }
        }
#endif
    }
    return 1;
}



void printStatsOnCriticalNodes(Node root, int *nbCriticalNodes) {
    makeStatsCriticalNodes(root, nbCriticalNodes);
//for (int i = 0; i < root->height; i++) printf("[%d] %d\n", i, nbCriticalNodes[i]);
    printf("[%d] %d\n", 0, nbCriticalNodes[0]);
    for (int i = 1; i < root->height; i++) if (nbCriticalNodes[i] != nbCriticalNodes[i-1]) printf("[%d] %d\n", i, nbCriticalNodes[i]);
    printf("[%d] %d\n", root->height-1, nbCriticalNodes[root->height-1]);
}