#include "treeordercomputationilp.h"
#include "unionfind.h"

void someLeaf(vector<size_t> &toLeaf, size_t nodeIndex, TemporalTree &tree)
{
    if (tree.isLeaf(nodeIndex))
    {
        toLeaf[nodeIndex] = nodeIndex;
        return;
    }
    size_t minChild = INT_MAX;
    for (auto child : tree.edgesHierarchy[nodeIndex])
    {
        someLeaf(toLeaf, child, tree);
        toLeaf[nodeIndex] = toLeaf[child];
    }
}

TemporalTree getILPSolution(GRBModel& model, TemporalTree tree) {
    TemporalTree outTree = tree.getSelfDeaggregated();
    outTree.computeReverseEdges();
    outTree.order.clear();
    map<int, vector<size_t>> t_nodes;
    map<int, vector<size_t>> t_leaves;
    vector<bool> isCurrentLeafDescendant(outTree.nodes.size(), false);
    for (size_t i = 1; i < outTree.nodes.size(); i++)
    {
        t_nodes[outTree.nodes[i].startTime()].push_back(i);
        if (outTree.isLeaf(i))
            t_leaves[outTree.nodes[i].startTime()].push_back(i);
    }
    for (auto const &[time, nodesattime] : t_nodes)
    {
        vector<size_t> leavesattime = t_leaves[time];
        vector<size_t> orderatt(leavesattime.size(), 0);
        for (int i = 0; i < leavesattime.size(); i++)
        {
            int cntbefore = 0;
            for (int j = 0; j < leavesattime.size(); j++)
            {
                if (i == j)
                    continue;
                size_t leaf1 = leavesattime[i];
                size_t leaf2 = leavesattime[j];
                stringstream ss1;
                stringstream ss2;
                ss1<<leaf1<<"<"<<leaf2;
                ss1<<leaf2<<"<"<<leaf1;
                

                if (leaf1 < leaf2 && model.getVarByName(ss1.str()).get(GRB_DoubleAttr_X) < 0.5)
                {
                    cntbefore++;
                }
                else if (leaf1 > leaf2 && model.getVarByName(ss2.str()).get(GRB_DoubleAttr_X) > 0.5)
                {
                    cntbefore++;
                }
            }
            orderatt[cntbefore] = leavesattime[i];
        }
        outTree.order.insert(outTree.order.end(), orderatt.begin(), orderatt.end());
        /*cout << "Order at time " << time << " is: ";
        for (auto nodeindex : orderatt)
        {
            cout << nodeindex << " ";
        }
        cout << endl;*/
    }
    return outTree;
}

TemporalTree computeTreeOrderILP(TemporalTree tree, GRBModel& model)
{
    auto starttime = high_resolution_clock::now();
    TemporalTree outTree = tree.getSelfDeaggregated();
    outTree.computeReverseEdges();


    vector<int> preOrderIndex(outTree.nodes.size());
    vector<int> postOrderIndex(outTree.nodes.size());
    int cnter = 0;
    treeorder::fillPreorder(preOrderIndex, cnter, outTree, 0);
    cnter = 0;
    treeorder::fillPostorder(postOrderIndex, cnter, outTree, 0);
    vector<size_t> toLeaf(outTree.nodes.size());
    someLeaf(toLeaf, 0, outTree);

    unordered_map<pair<size_t, size_t>, GRBVar, pairhash> ordvars;

    map<int, vector<size_t>> t_nodes;
    map<int, vector<size_t>> t_leaves;
    vector<bool> isCurrentLeafDescendant(outTree.nodes.size(), false);
    for (size_t i = 1; i < outTree.nodes.size(); i++)
    {
        t_nodes[outTree.nodes[i].startTime()].push_back(i);
        if (outTree.isLeaf(i))
            t_leaves[outTree.nodes[i].startTime()].push_back(i);
    }

    for (auto const &[time, nodesattime] : t_nodes)
    {
        // cout<<nodesattime.size()<<endl;
        vector<size_t> leavesattime = t_leaves[time];

        // create ordering variables
        for (int i = 0; i < leavesattime.size(); i++)
        {
            for (int j = i + 1; j < leavesattime.size(); j++)
            {
                stringstream ss;
                ss<<leavesattime[i]<<"<"<<leavesattime[j];
                ordvars[{leavesattime[i], leavesattime[j]}] = model.addVar(0, 1, 0, GRB_BINARY, ss.str());
            }
        }

        // transitivity constraints
        for (int i = 0; i < leavesattime.size(); i++)
        {
            for (int j = i + 1; j < leavesattime.size(); j++)
            {
                for (int k = j + 1; k < leavesattime.size(); k++)
                {
                    auto v1 = ordvars[{leavesattime[i], leavesattime[j]}];
                    auto v2 = ordvars[{leavesattime[j], leavesattime[k]}];
                    auto v3 = ordvars[{leavesattime[i], leavesattime[k]}];
                    model.addConstr(0 <= v1 + v2 - v3);
                    model.addConstr(v1 + v2 - v3 <= 1);
                }
            }
        }

        // tree constraints
        for (const auto &nodeIndex : nodesattime)
        {
            if (!outTree.isLeaf(nodeIndex) && outTree.edgesHierarchy[nodeIndex].size() > 1)
            {
                vector<size_t> leaves;
                outTree.computeLeafDescendants(nodeIndex, leaves);
                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = true;

                for (const auto &leafoutside : leavesattime)
                {
                    if (isCurrentLeafDescendant[leafoutside])
                        continue;

                    // now force the leaf to be really outside of the current subtree
                    for (int i = 0; i < leaves.size(); i++)
                    {
                        for (int j = i + 1; j < leaves.size(); j++)
                        {
                            auto leafinside1 = min(leaves[i], leaves[j]);
                            auto leafinside2 = max(leaves[i], leaves[j]);

                            if (leafoutside < leafinside1 && leafoutside < leafinside2)
                            {
                                model.addConstr(ordvars[{leafoutside, leafinside1}] == ordvars[{leafoutside, leafinside2}]);
                            }
                            else if (leafoutside > leafinside1 && leafoutside > leafinside2)
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == ordvars[{leafinside2, leafoutside}]);
                            }
                            else
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == 1 - ordvars[{leafoutside, leafinside2}]);
                            }
                        }
                    }
                }

                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = false;
            }
        }
    }

    // crossing constraints
    for (auto const &[time, nodesattime] : t_nodes)
    {
        // edges that go from one time to the next
        vector<pair<size_t, size_t>> edges;
        for (auto v : nodesattime)
        {
            if (outTree.edgesTime.find(v) != outTree.edgesTime.end())
            {
                for (auto w : outTree.edgesTime[v])
                {
                    edges.push_back({v, w});
                }
            }
        }
        // now check for every pair of edges if they cross
        for (size_t i = 0; i < edges.size(); i++)
        {
            for (size_t j = i + 1; j < edges.size(); j++)
            {
                auto e1 = edges[i];
                auto e2 = edges[j];

                if (e1.first == e2.first || e1.second == e2.second)
                    continue;
                // we dont want do count those crossings
                if (treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.first, e2.first) || treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.second, e2.second))
                    continue;

                auto leaf1left = toLeaf[e1.first];
                auto leaf1right = toLeaf[e1.second];
                auto leaf2left = toLeaf[e2.first];
                auto leaf2right = toLeaf[e2.second];

                pair<size_t, size_t> varleftind = {min(leaf1left, leaf2left), max(leaf1left, leaf2left)};
                pair<size_t, size_t> varrightind = {min(leaf1right, leaf2right), max(leaf1right, leaf2right)};

                assert(ordvars.find(varleftind) != ordvars.end());
                assert(ordvars.find(varrightind) != ordvars.end());

                auto varleft = ordvars[varleftind];
                auto varright = ordvars[varrightind];

                auto varCrossing = model.addVar(0, 1, 1, GRB_BINARY);
                // now set crossing variable
                if ((leaf1left < leaf2left) == (leaf1right < leaf2right))
                {
                    // if the two variables have different values then we have a crossing
                    model.addConstr(varCrossing >= varleft - varright);
                    model.addConstr(varCrossing >= varright - varleft);
                }
                else
                {
                    model.addConstr(varCrossing >= 1 - varleft - varright);
                    model.addConstr(varCrossing >= varright - 1 + varleft);
                }
            }
        }
    }
    model.set(GRB_IntAttr_ModelSense, GRB_MINIMIZE);
    auto t2 = high_resolution_clock::now();
    duration<double, std::milli> ms_double = t2 - starttime;
    model.getEnv().set(GRB_DoubleParam_TimeLimit, (double)TIMELIMIT.count() - ms_double.count()/((double)1000));
    model.optimize();

    if (model.get(GRB_IntAttr_SolCount) > 0)
    {
        outTree.order.clear();
        for (auto const &[time, nodesattime] : t_nodes)
        {
            vector<size_t> leavesattime = t_leaves[time];
            vector<size_t> orderatt(leavesattime.size(), 0);
            for (int i = 0; i < leavesattime.size(); i++)
            {
                int cntbefore = 0;
                for (int j = 0; j < leavesattime.size(); j++)
                {
                    if (i == j)
                        continue;
                    size_t leaf1 = leavesattime[i];
                    size_t leaf2 = leavesattime[j];

                    if (leaf1 < leaf2 && ordvars[{leaf1, leaf2}].get(GRB_DoubleAttr_X) < 0.5)
                    {
                        cntbefore++;
                    }
                    else if (leaf1 > leaf2 && ordvars[{leaf2, leaf1}].get(GRB_DoubleAttr_X) > 0.5)
                    {
                        cntbefore++;
                    }
                }
                orderatt[cntbefore] = leavesattime[i];
            }
            outTree.order.insert(outTree.order.end(), orderatt.begin(), orderatt.end());
            /*cout << "Order at time " << time << " is: ";
            for (auto nodeindex : orderatt)
            {
                cout << nodeindex << " ";
            }
            cout << endl;*/
        }
    }
    else
    {
        throw GRBException("No solution found");
    }
    return outTree;
}

TemporalTree computeTreeOrderILPQuadr(TemporalTree tree, GRBModel& model) 
{
    auto starttime = high_resolution_clock::now();
    TemporalTree outTree = tree.getSelfDeaggregated();
    outTree.computeReverseEdges();

    vector<int> preOrderIndex(outTree.nodes.size());
    vector<int> postOrderIndex(outTree.nodes.size());
    int cnter = 0;
    treeorder::fillPreorder(preOrderIndex, cnter, outTree, 0);
    cnter = 0;
    treeorder::fillPostorder(postOrderIndex, cnter, outTree, 0);
    vector<size_t> toLeaf(outTree.nodes.size());
    someLeaf(toLeaf, 0, outTree);

    unordered_map<pair<size_t, size_t>, GRBVar, pairhash> ordvars;

    map<int, vector<size_t>> t_nodes;
    map<int, vector<size_t>> t_leaves;
    vector<bool> isCurrentLeafDescendant(outTree.nodes.size(), false);
    for (size_t i = 1; i < outTree.nodes.size(); i++)
    {
        t_nodes[outTree.nodes[i].startTime()].push_back(i);
        if (outTree.isLeaf(i))
            t_leaves[outTree.nodes[i].startTime()].push_back(i);
    }

    for (auto const &[time, nodesattime] : t_nodes)
    {
        // cout<<nodesattime.size()<<endl;
        vector<size_t> leavesattime = t_leaves[time];

        // create ordering variables
        for (int i = 0; i < leavesattime.size(); i++)
        {
            for (int j = i + 1; j < leavesattime.size(); j++)
            {
                stringstream ss;
                ss<<leavesattime[i]<<"<"<<leavesattime[j];
                ordvars[{leavesattime[i], leavesattime[j]}] = model.addVar(0, 1, 0, GRB_BINARY, ss.str());
            }
        }

        // transitivity constraints
        for (int i = 0; i < leavesattime.size(); i++)
        {
            for (int j = i + 1; j < leavesattime.size(); j++)
            {
                for (int k = j + 1; k < leavesattime.size(); k++)
                {
                    auto v1 = ordvars[{leavesattime[i], leavesattime[j]}];
                    auto v2 = ordvars[{leavesattime[j], leavesattime[k]}];
                    auto v3 = ordvars[{leavesattime[i], leavesattime[k]}];
                    model.addConstr(0 <= v1 + v2 - v3);
                    model.addConstr(v1 + v2 - v3 <= 1);
                }
            }
        }

        // tree constraints
        for (const auto &nodeIndex : nodesattime)
        {
            if (!outTree.isLeaf(nodeIndex) && outTree.edgesHierarchy[nodeIndex].size() > 1)
            {
                vector<size_t> leaves;
                outTree.computeLeafDescendants(nodeIndex, leaves);
                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = true;

                for (const auto &leafoutside : leavesattime)
                {
                    if (isCurrentLeafDescendant[leafoutside])
                        continue;

                    // now force the leaf to be really outside of the current subtree
                    for (int i = 0; i < leaves.size(); i++)
                    {
                        for (int j = i + 1; j < leaves.size(); j++)
                        {
                            auto leafinside1 = min(leaves[i], leaves[j]);
                            auto leafinside2 = max(leaves[i], leaves[j]);

                            if (leafoutside < leafinside1 && leafoutside < leafinside2)
                            {
                                model.addConstr(ordvars[{leafoutside, leafinside1}] == ordvars[{leafoutside, leafinside2}]);
                            }
                            else if (leafoutside > leafinside1 && leafoutside > leafinside2)
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == ordvars[{leafinside2, leafoutside}]);
                            }
                            else
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == 1 - ordvars[{leafoutside, leafinside2}]);
                            }
                        }
                    }
                }

                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = false;
            }
        }
    }

    // crossing constraints
    GRBQuadExpr obj = 0.0;
    for (auto const &[time, nodesattime] : t_nodes)
    {
        // edges that go from one time to the next
        vector<pair<size_t, size_t>> edges;
        for (auto v : nodesattime)
        {
            if (outTree.edgesTime.find(v) != outTree.edgesTime.end())
            {
                for (auto w : outTree.edgesTime[v])
                {
                    edges.push_back({v, w});
                }
            }
        }
        // now check for every pair of edges if they cross
        for (size_t i = 0; i < edges.size(); i++)
        {
            for (size_t j = i + 1; j < edges.size(); j++)
            {
                auto e1 = edges[i];
                auto e2 = edges[j];

                if (e1.first == e2.first || e1.second == e2.second)
                    continue;
                // we dont want do count those crossings
                if (treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.first, e2.first) || treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.second, e2.second))
                    continue;

                auto leaf1left = toLeaf[e1.first];
                auto leaf1right = toLeaf[e1.second];
                auto leaf2left = toLeaf[e2.first];
                auto leaf2right = toLeaf[e2.second];

                pair<size_t, size_t> varleftind = {min(leaf1left, leaf2left), max(leaf1left, leaf2left)};
                pair<size_t, size_t> varrightind = {min(leaf1right, leaf2right), max(leaf1right, leaf2right)};

                assert(ordvars.find(varleftind) != ordvars.end());
                assert(ordvars.find(varrightind) != ordvars.end());

                auto varleft = ordvars[varleftind];
                auto varright = ordvars[varrightind];

                auto varCrossing = model.addVar(0, 1, 1, GRB_BINARY);
                // now set crossing variable
                if ((leaf1left < leaf2left) == (leaf1right < leaf2right))
                {
                    // if the two variables have different values then we have a crossing
                    obj += varleft*(1-varright);
                    obj += (1-varleft)*varright;
                }
                else
                {
                    obj += varleft*varright;
                    obj += (1-varleft)*(1-varright);
                }
            }
        }
    }
    model.set(GRB_IntAttr_ModelSense, GRB_MINIMIZE);
    auto t2 = high_resolution_clock::now();
    duration<double, std::milli> ms_double = t2 - starttime;
    model.getEnv().set(GRB_DoubleParam_TimeLimit, (double)TIMELIMIT.count() - ms_double.count()/((double)1000));
    model.optimize();

    if (model.get(GRB_IntAttr_SolCount) > 0)
    {
        outTree.order.clear();
        for (auto const &[time, nodesattime] : t_nodes)
        {
            vector<size_t> leavesattime = t_leaves[time];
            vector<size_t> orderatt(leavesattime.size(), 0);
            for (int i = 0; i < leavesattime.size(); i++)
            {
                int cntbefore = 0;
                for (int j = 0; j < leavesattime.size(); j++)
                {
                    if (i == j)
                        continue;
                    size_t leaf1 = leavesattime[i];
                    size_t leaf2 = leavesattime[j];

                    if (leaf1 < leaf2 && ordvars[{leaf1, leaf2}].get(GRB_DoubleAttr_X) < 0.5)
                    {
                        cntbefore++;
                    }
                    else if (leaf1 > leaf2 && ordvars[{leaf2, leaf1}].get(GRB_DoubleAttr_X) > 0.5)
                    {
                        cntbefore++;
                    }
                }
                orderatt[cntbefore] = leavesattime[i];
            }
            outTree.order.insert(outTree.order.end(), orderatt.begin(), orderatt.end());
            /*cout << "Order at time " << time << " is: ";
            for (auto nodeindex : orderatt)
            {
                cout << nodeindex << " ";
            }
            cout << endl;*/
        }
    }
    else
    {
        throw GRBException("No solution found");
    }
    return outTree;
}

TemporalTree computeTreeOrderLeafEdgesILP(TemporalTree tree, GRBModel& model)
{
    auto starttime = high_resolution_clock::now();
    TemporalTree outTree = tree.getSelfDeaggregated();
    outTree.computeReverseEdges();

    vector<int> preOrderIndex(outTree.nodes.size());
    vector<int> postOrderIndex(outTree.nodes.size());
    int cnter = 0;
    treeorder::fillPreorder(preOrderIndex, cnter, outTree, 0);
    cnter = 0;
    treeorder::fillPostorder(postOrderIndex, cnter, outTree, 0);
    vector<size_t> toLeaf(outTree.nodes.size());
    someLeaf(toLeaf, 0, outTree);

    unordered_map<pair<size_t, size_t>, GRBVar, pairhash> ordvars;

    map<int, vector<size_t>> t_nodes;
    map<int, vector<size_t>> t_leaves;
    vector<bool> isCurrentLeafDescendant(outTree.nodes.size(), false);
    for (size_t i = 1; i < outTree.nodes.size(); i++)
    {
        t_nodes[outTree.nodes[i].startTime()].push_back(i);
        if (outTree.isLeaf(i))
            t_leaves[outTree.nodes[i].startTime()].push_back(i);
    }

    for (auto const &[time, nodesattime] : t_nodes)
    {
        // cout<<nodesattime.size()<<endl;
        vector<size_t> leavesattime = t_leaves[time];

        // create ordering variables
        for (int i = 0; i < leavesattime.size(); i++)
        {
            for (int j = i + 1; j < leavesattime.size(); j++)
            {
                stringstream ss;
                ss<<leavesattime[i]<<"<"<<leavesattime[j];
                ordvars[{leavesattime[i], leavesattime[j]}] = model.addVar(0, 1, 0, GRB_BINARY, ss.str());
            }
        }

        // transitivity constraints
        for (int i = 0; i < leavesattime.size(); i++)
        {
            for (int j = i + 1; j < leavesattime.size(); j++)
            {
                for (int k = j + 1; k < leavesattime.size(); k++)
                {
                    auto v1 = ordvars[{leavesattime[i], leavesattime[j]}];
                    auto v2 = ordvars[{leavesattime[j], leavesattime[k]}];
                    auto v3 = ordvars[{leavesattime[i], leavesattime[k]}];
                    model.addConstr(0 <= v1 + v2 - v3);
                    model.addConstr(v1 + v2 - v3 <= 1);
                }
            }
        }

        // tree constraints
        for (const auto &nodeIndex : nodesattime)
        {
            if (!outTree.isLeaf(nodeIndex) && outTree.edgesHierarchy[nodeIndex].size() > 1)
            {
                vector<size_t> leaves;
                outTree.computeLeafDescendants(nodeIndex, leaves);
                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = true;

                for (const auto &leafoutside : leavesattime)
                {
                    if (isCurrentLeafDescendant[leafoutside])
                        continue;

                    // now force the leaf to be really outside of the current subtree
                    for (int i = 0; i < leaves.size(); i++)
                    {
                        for (int j = i + 1; j < leaves.size(); j++)
                        {
                            auto leafinside1 = min(leaves[i], leaves[j]);
                            auto leafinside2 = max(leaves[i], leaves[j]);

                            if (leafoutside < leafinside1 && leafoutside < leafinside2)
                            {
                                model.addConstr(ordvars[{leafoutside, leafinside1}] == ordvars[{leafoutside, leafinside2}]);
                            }
                            else if (leafoutside > leafinside1 && leafoutside > leafinside2)
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == ordvars[{leafinside2, leafoutside}]);
                            }
                            else
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == 1 - ordvars[{leafoutside, leafinside2}]);
                            }
                        }
                    }
                }

                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = false;
            }
        }
    }

    // crossing constraints

    outTree.computeLeafEdges();
    for (auto const &[time, nodesattime] : t_nodes)
    {
        // edges that go from one time to the next
        vector<pair<size_t, size_t>> edges;
        for (auto v : nodesattime)
        {
            if (outTree.leafEdges.find(v) != outTree.leafEdges.end())
            {
                for (auto w : outTree.leafEdges[v])
                {
                    edges.push_back({v, w});
                }
            }
        }
        // now check for every pair of edges if they cross
        for (size_t i = 0; i < edges.size(); i++)
        {
            for (size_t j = i + 1; j < edges.size(); j++)
            {
                auto e1 = edges[i];
                auto e2 = edges[j];

                if (e1.first == e2.first || e1.second == e2.second)
                    continue;
                // we dont want do count those crossings
                if (treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.first, e2.first) || treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.second, e2.second))
                    continue;

                auto leaf1left = toLeaf[e1.first];
                auto leaf1right = toLeaf[e1.second];
                auto leaf2left = toLeaf[e2.first];
                auto leaf2right = toLeaf[e2.second];

                pair<size_t, size_t> varleftind = {min(leaf1left, leaf2left), max(leaf1left, leaf2left)};
                pair<size_t, size_t> varrightind = {min(leaf1right, leaf2right), max(leaf1right, leaf2right)};

                assert(ordvars.find(varleftind) != ordvars.end());
                assert(ordvars.find(varrightind) != ordvars.end());

                auto varleft = ordvars[varleftind];
                auto varright = ordvars[varrightind];

                auto varCrossing = model.addVar(0, 1, 1, GRB_BINARY);
                // now set crossing variable
                if ((leaf1left < leaf2left) == (leaf1right < leaf2right))
                {
                    // if the two variables have different values then we have a crossing
                    model.addConstr(varCrossing >= varleft - varright);
                    model.addConstr(varCrossing >= varright - varleft);
                }
                else
                {
                    model.addConstr(varCrossing >= 1 - varleft - varright);
                    model.addConstr(varCrossing >= varright - 1 + varleft);
                }
            }
        }
    }
    model.set(GRB_IntAttr_ModelSense, GRB_MINIMIZE);
    auto t2 = high_resolution_clock::now();
    duration<double, std::milli> ms_double = t2 - starttime;
    model.getEnv().set(GRB_DoubleParam_TimeLimit, (double)TIMELIMIT.count() - ms_double.count()/((double)1000));
    model.optimize();

    if (model.get(GRB_IntAttr_SolCount) > 0)
    {
        outTree.order.clear();
        for (auto const &[time, nodesattime] : t_nodes)
        {
            vector<size_t> leavesattime = t_leaves[time];
            vector<size_t> orderatt(leavesattime.size(), 0);
            for (int i = 0; i < leavesattime.size(); i++)
            {
                int cntbefore = 0;
                for (int j = 0; j < leavesattime.size(); j++)
                {
                    if (i == j)
                        continue;
                    size_t leaf1 = leavesattime[i];
                    size_t leaf2 = leavesattime[j];

                    if (leaf1 < leaf2 && ordvars[{leaf1, leaf2}].get(GRB_DoubleAttr_X) < 0.5)
                    {
                        cntbefore++;
                    }
                    else if (leaf1 > leaf2 && ordvars[{leaf2, leaf1}].get(GRB_DoubleAttr_X) > 0.5)
                    {
                        cntbefore++;
                    }
                }
                orderatt[cntbefore] = leavesattime[i];
            }
            outTree.order.insert(outTree.order.end(), orderatt.begin(), orderatt.end());
            /*cout << "Order at time " << time << " is: ";
            for (auto nodeindex : orderatt)
            {
                cout << nodeindex << " ";
            }
            cout << endl;*/
        }
    }
    else
    {
        throw GRBException("No solution found");
    }
    return outTree;
}

/*
void recursive_aggregation(TemporalTree& tree, size_t nodeIndex, vector<AggrType>& aggrtypes, union_find& DSU, unordered_map<int,int>& to_ind_union) {
    if (aggrtypes[to_ind_union[nodeIndex]] != AggrType::NOT_SET) return;
    bool isChildrenAggregated = true;
    if (!tree.isLeaf(nodeIndex)) {
        for (auto child: tree.edgesHierarchy[nodeIndex]) {
            recursive_aggregation(tree, child, aggrtypes, DSU, to_ind_union);
            if (aggrtypes[to_ind_union[child]] != AggrType::YES) isChildrenAggregated = false;
        }
    }

    if (isChildrenAggregated) {
        // check if I am also aggregated, if yes, union_find all children and me
        if (tree.reverseEdgesTime.find(nodeIndex) != tree.reverseEdgesTime.end() && tree.reverseEdgesTime[nodeIndex].size() == 1 &&
            tree.edgesTime[tree.reverseEdgesTime[nodeIndex][0]].size() == 1) {
            // exactly one time edge from me from predecessor
            if (!tree.isLeaf(nodeIndex)) {
                // merge with all children
                for (auto child: tree.edgesHierarchy[nodeIndex]) {
                    DSU.unite(to_ind_union[nodeIndex], to_ind_union[child]);
                }
            }
            aggrtypes[to_ind_union[nodeIndex]] = AggrType::YES;
        }
        else {
            aggrtypes[to_ind_union[nodeIndex]] = AggrType::NO;
        }
    }
    else {
        aggrtypes[to_ind_union[nodeIndex]] = AggrType::NO;
    }
}

TemporalTree computeTreeOrderILPAggregated(TemporalTree tree)
{
    TemporalTree outTree = tree.getSelfDeaggregated();
    outTree.checkConsistency();
    outTree.computeReverseEdges();

    vector<int> preOrderIndex(outTree.nodes.size());
    vector<int> postOrderIndex(outTree.nodes.size());
    int cnter = 0;
    treeorder::fillPreorder(preOrderIndex, cnter, outTree, 0);
    cnter = 0;
    treeorder::fillPostorder(postOrderIndex, cnter, outTree, 0);
    vector<size_t> toLeaf(outTree.nodes.size());
    someLeaf(toLeaf, 0, outTree);

    GRBEnv env(true);
    // env.set(GRB_IntParam_Threads, 1);
    // env.set(GRB_IntParam_OutputFlag, 0);
    env.set(GRB_DoubleParam_MemLimit, 4);
    env.start();
    GRBModel model(env);

    unordered_map<pair<size_t, size_t>, GRBVar, pairhash> ordvars;

    map<int, vector<size_t>> t_nodes;
    map<int, vector<size_t>> t_leaves;
    unordered_map<int, union_find> t_unions;
    unordered_map<int, unordered_map<int,int>> to_ind_union;
    vector<bool> isCurrentLeafDescendant(outTree.nodes.size(), false);
    for (size_t i = 1; i < outTree.nodes.size(); i++)
    {
        t_nodes[outTree.nodes[i].startTime()].push_back(i);
        if (outTree.isLeaf(i))
            t_leaves[outTree.nodes[i].startTime()].push_back(i);
    }

    for (auto const &[time, nodesattime] : t_nodes)
    {
        t_unions[time] = union_find(nodesattime.size());
        for (int i = 0; i < nodesattime.size(); i++) to_ind_union[time][nodesattime[i]] = i;
        vector<size_t> leavesattime = t_leaves[time];
        vector<AggrType> aggrtypes(nodesattime.size(), AggrType::NOT_SET);
        for (auto nodeIndex: nodesattime) recursive_aggregation(outTree, nodeIndex, aggrtypes, t_unions[time], to_ind_union[time]);
        unordered_set<int> representatives;
        vector<size_t> leafs;
        for (auto nodeIndex: nodesattime) {
            if (outTree.isLeaf(nodeIndex)) {
                representatives.insert(t_unions[time].find(to_ind_union[time][nodeIndex]));
                leafs.push_back(nodeIndex);
            }
        }

        vector<int> representativesvec(representatives.begin(), representatives.end());
        cout<<representatives.size()<<" "<<leafs.size()<<endl;
        // create ordering variables
        for (int i = 0; i < representativesvec.size(); i++)
        {
            for (int j = i + 1; j < representativesvec.size(); j++)
            {
                ordvars[{representativesvec[i], representativesvec[j]}] = model.addVar(0, 1, 0, GRB_BINARY);
            }
        }

        // transitivity constraints
        for (int i = 0; i < representativesvec.size(); i++)
        {
            for (int j = i + 1; j < representativesvec.size(); j++)
            {
                for (int k = j + 1; k < representativesvec.size(); k++)
                {
                    auto v1 = ordvars[{representativesvec[i], representativesvec[j]}];
                    auto v2 = ordvars[{representativesvec[j], representativesvec[k]}];
                    auto v3 = ordvars[{representativesvec[i], representativesvec[k]}];
                    model.addConstr(0 <= v1 + v2 - v3);
                    model.addConstr(v1 + v2 - v3 <= 1);
                }
            }
        }

        // tree constraints
        /*for (const auto &nodeIndex : nodesattime)
        {
            if (!outTree.isLeaf(nodeIndex) && outTree.edgesHierarchy[nodeIndex].size() > 1)
            {
                vector<size_t> leaves;
                computeLeafDescendants(nodeIndex, outTree, leaves);
                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = true;

                for (const auto &leafoutside : leavesattime)
                {
                    if (isCurrentLeafDescendant[leafoutside])
                        continue;

                    // now force the leaf to be really outside of the current subtree
                    for (int i = 0; i < leaves.size(); i++)
                    {
                        for (int j = i + 1; j < leaves.size(); j++)
                        {
                            auto leafinside1 = min(leaves[i], leaves[j]);
                            auto leafinside2 = max(leaves[i], leaves[j]);

                            if (leafoutside < leafinside1 && leafoutside < leafinside2)
                            {
                                model.addConstr(ordvars[{leafoutside, leafinside1}] == ordvars[{leafoutside, leafinside2}]);
                            }
                            else if (leafoutside > leafinside1 && leafoutside > leafinside2)
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == ordvars[{leafinside2, leafoutside}]);
                            }
                            else
                            {
                                model.addConstr(ordvars[{leafinside1, leafoutside}] == 1 - ordvars[{leafoutside, leafinside2}]);
                            }
                        }
                    }
                }

                for (const auto &leaf : leaves)
                    isCurrentLeafDescendant[leaf] = false;
            }
        }
    }

    // crossing constraints
    for (auto const &[time, nodesattime] : t_nodes)
    {
        // edges that go from one time to the next
        vector<pair<size_t, size_t>> edges;
        for (auto v : nodesattime)
        {
            for (auto w : outTree.edgesTime[v])
            {
                edges.push_back({v, w});
            }
        }
        // now check for every pair of edges if they cross
        for (size_t i = 0; i < edges.size(); i++)
        {
            for (size_t j = i + 1; j < edges.size(); j++)
            {
                auto e1 = edges[i];
                auto e2 = edges[j];

                if (e1.first == e2.first || e1.second == e2.second)
                    continue;
                // we dont want do count those crossings
                if (treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.first, e2.first) || treeorder::isChildDescRelation(preOrderIndex, postOrderIndex, e1.second, e2.second))
                    continue;

                auto leaf1left = toLeaf[e1.first];
                auto leaf1right = toLeaf[e1.second];
                auto leaf2left = toLeaf[e2.first];
                auto leaf2right = toLeaf[e2.second];

                auto varleft = ordvars[{min(leaf1left, leaf2left), max(leaf1left, leaf2left)}];
                auto varright = ordvars[{min(leaf1right, leaf2right), max(leaf1right, leaf2right)}];

                auto varCrossing = model.addVar(0, 1, 1, GRB_BINARY);
                // now set crossing variable
                if ((leaf1left < leaf2left) == (leaf1right < leaf2right))
                {
                    // if the two variables have different values then we have a crossing
                    model.addConstr(varCrossing >= varleft-varright);
                    model.addConstr(varCrossing >= varright-varleft);
                }
                else
                {
                    model.addConstr(varCrossing >= 1-varleft-varright);
                    model.addConstr(varCrossing >= varright-1+varleft);
                }
            }
        }
    }
    model.set(GRB_IntAttr_ModelSense, GRB_MINIMIZE);
    model.optimize();

    if (model.get(GRB_IntAttr_SolCount) > 0)
    {
        cout<<"is optimal"<<endl;
        outTree.order.clear();
        for (auto const &[time, nodesattime] : t_nodes)
        {
            vector<size_t> leavesattime = t_leaves[time];
            vector<size_t> orderatt(leavesattime.size(), 0);
            for (int i = 0; i < leavesattime.size(); i++)
            {
                int cntbefore = 0;
                for (int j = 0; j < leavesattime.size(); j++)
                {
                    if (i == j)
                        continue;
                    size_t leaf1 = leavesattime[i];
                    size_t leaf2 = leavesattime[j];

                    if (leaf1 < leaf2 && ordvars[{leaf1, leaf2}].get(GRB_DoubleAttr_X) < 0.5)
                    {
                        cntbefore++;
                    }
                    else if (leaf1 > leaf2 && ordvars[{leaf2, leaf1}].get(GRB_DoubleAttr_X) > 0.5)
                    {
                        cntbefore++;
                    }
                }
                orderatt[cntbefore] = leavesattime[i];
            }
            outTree.order.insert(outTree.order.end(), orderatt.begin(), orderatt.end());
            cout << "Order at time " << time << " is: ";
            for (auto nodeindex : orderatt)
            {
                cout << outTree.nodes[nodeindex].name << " ";
            }
            cout << endl;
        }
        assert(treeorder::fitsWithTree(outTree, outTree.order));
        assert(getCntViolatedHierarchyConstraints(outTree) == 0);
        cout << "Number of Crossings: " << treeorder::countCrossings(outTree) << endl;
        // outTree.drawTree();
    }
    else
    {
        throw GRBException("No solution found");
    }

    return tree;
}*/