discretizationsolver.cpp

/********************************************************************************************
 Name:       DiscretizationSolver
 Discrete optimization for graph discretization
 Author:     J.Omer
 Sources:    C++
 License:    GNU General Public License v.2
 History:
 *********************************************************************************************/

#include "discretizationsolver.hpp"
#include <ctime>



//

ILOLAZYCONSTRAINTCALLBACK2(CyclesLazyConstraints, Instance &, inst, DiscretizationSolver&, solver) {

#ifdef VERBOSE
    std::cout << "revorder: check lazy dicycle constraints" << std::endl;
#endif

    // Search for cycles in the current solution digraph
    //
    // generate the adjacency lists of the vertices corresponding to the current
    // cplex solution
    float tolerance = 1e-06;
    std::map<Vertex*, std::vector<Vertex*> > adjlist;
    std::map<Vertex*, int > nbrefs;
    for (Vertex* v : inst.vertices_) {
        adjlist[v] = std::vector<Vertex*>();
        nbrefs[v] = 0;
    }
    for (Vertex* u : inst.vertices_) {
        if (solver.modeltype_ == WITNESS) {
            if (getValue(solver.isvertexinclique(u) >= 1 - tolerance)) continue;
        }
        for (Vertex* v : u->neighbors_) {
            if (getValue(solver.isbefore(u, v)) >= 1 - tolerance) {
                adjlist[u].push_back(v);
                nbrefs[v]++;
            }
        }
    }
    //
    // search for the root of the digraph
    Vertex* root = nullptr;
    int minnbrefs = inst.nbvertices_;
    if (solver.modeltype_ == WITNESS) {
        for (Vertex* v : inst.vertices_) {
            if (getValue(solver.isvertexinclique(v) >= 1 - tolerance)) {
                root = v;
                break;
            }
        }
    } else {
        for (Vertex* v : inst.vertices_) {
            if (nbrefs[v] < minnbrefs) {
                minnbrefs = nbrefs[v];
                root = v;
                if (minnbrefs < inst.L()) break;
            }
        }
    }
    //
    // call the enumeration of cycles, returns false if there are none
    std::vector<std::vector<Vertex*> > cycles;
    bool iscyclic = solver.enumeratecycles(adjlist, root, cycles);

    // Look for dicycle inequalities : we restrict the search to the dicycles
    // that contain at least one edge in the distance graph
    // besides that, the search is performed using brute force
    int nbaddedcuts = 0;
    IloEnv env = getEnv();
    if (iscyclic) {
        for (std::vector<Vertex*> path : cycles) {
            IloExpr sumedges(env);
            int cyclelength = path.size();
            Vertex* u = path.back();
            for (int i = 0; i < cyclelength; i++) {
                Vertex* v = path[i];
                sumedges += solver.isbefore(u, v);
                u = v;
            }
            if ((solver.modeltype_ == WITNESS) && (cyclelength <= inst.L() + 1)) {
                add(sumedges - solver.isvertexinclique(path.front()) <= cyclelength - 1).end();
                nbaddedcuts++;
            } else {
                add(sumedges <= cyclelength - 1).end();
                nbaddedcuts++;
            }
        }
#ifdef VERBOSE
        std::cout << "revorder: number of violated dicycle constraints : " << nbaddedcuts << std::endl;
#endif
    } else {
#ifdef VERBOSE
        std::cout << "revorder: the graph is acyclic " << std::endl;
#endif
    }
    return;
}

// User branch callback that gives priority to branching on clique variables

ILOBRANCHCALLBACK2(BranchOnCliqueVariables, DiscretizationSolver&, solver, IloCplex::MIPCallbackI::NodeId&, nodeId) {
    if (getBranchType() != BranchOnVariable)
        return;

    if (nodeId == getNodeId())
        return;
    nodeId = getNodeId();

    // Branch on the fractionnary clique variable closest to 1
    IntegerFeasibilityArray feas;
    IloNumArray val;
    try {
        val = IloNumArray(getEnv());
        feas = IntegerFeasibilityArray(getEnv());
        getValues(val, solver.isclique_);
        getFeasibilities(feas, solver.isclique_);
        //
        // get the clique variable with maximum integer infeasibility
        IloInt bestvar = -1;
        IloNum maxval = 0.3;
        IloInt cols = solver.isclique_.getSize();
        for (IloInt j = 0; j < cols; j++) {
            if (feas[j] == Infeasible) {
                if (val[j] >= 1 - 1.0e-6) {
                    bestvar = -1;
                    break;
                }
                if (std::min(val[j], 1 - val[j]) > maxval) {
                    bestvar = j;
                    maxval = std::min(val[j], 1 - val[j]);
                }
            } else if (val[j] >= 1 - 1.0e-6) {
                bestvar = -1;
                break;
            }
        }
        //
        // if there is at least one fractionnary clique variable create two branches
        if (bestvar >= 0) {
            //            HeuristicNode* data = new HeuristicNode();
            //            data->fixcliqueindex(bestvar);
            makeBranch(solver.isclique_[bestvar], val[bestvar], IloCplex::BranchUp, getObjValue());
            makeBranch(solver.isclique_[bestvar], val[bestvar], IloCplex::BranchDown, getObjValue());
#ifdef VERBOSE
            std::cout << "revorder: clique " << bestvar << " bounds " << getLB(solver.isclique_[bestvar]) << " " << getUB(solver.isclique_[bestvar]) << ", value " << val[bestvar] << ", feasibility " << feas[bestvar] << std::endl;
            std::cout << "revorder: branch callback: branch on the variable isclique_" << bestvar << std::endl;
#endif
        }
        //        else {
        //            Vertex* bestu = NULL;
        //            Vertex* bestv = NULL;
        //            for (Vertex* u : solver.inst_->vertices_) {
        //                std::map<Vertex*, IloNum> isbeforeval;
        //                float sumrefs = 0;
        //                for (Vertex* v : u->neighbors_) {
        //                    isbeforeval[v] = getValue(solver.isbefore_[v][u]);
        //                    sumrefs += isbeforeval[v];
        //                }
        //                if (sumrefs == solver.inst_->U()) {
        //                    for (Vertex* v : u->neighbors_) {
        //                        if (getFeasibility(solver.isbefore_[v][u]) == Infeasible) {
        //                            if (std::min(isbeforeval[v], 1 - isbeforeval[v]) > maxval) {
        //                                bestu = u;
        //                                bestv = v;
        //                                maxval = std::min(isbeforeval[v], 1 - isbeforeval[v]);
        //                            }
        //                        }
        //                    }
        //                }
        //            }
        //            if (bestu) {
        //                // HeuristicNode* data = new HeuristicNode();
        //                // data->fixcliqueindex(bestvar);
        //                std::cout << "revorder: branch on isbefore_" << bestv->id_ << "_" << bestu->id_ << std::endl;
        //                makeBranch(solver.isbefore_[bestv][bestu], getValue(solver.isbefore_[bestv][bestu]), IloCplex::BranchUp, getObjValue());
        //                makeBranch(solver.isbefore_[bestv][bestu], getValue(solver.isbefore_[bestv][bestu]), IloCplex::BranchDown, getObjValue());
        //            }
        //        }
    } catch (...) {

        val.end();
        feas.end();
        throw;
    }
    val.end();
    feas.end();
}


//
// constructors and destructor

Clique::Clique(std::vector<Vertex *> vert) : nbvertices_(vert.size()), vertices_(vert) {
}

Clique::~Clique() {
}

DiscretizationSolver::DiscretizationSolver(Instance* inst, Algorithm algo) : algo_(algo), inst_(inst) {
    this->bestnbfullref_ = -1;
    this->bestcliqueid_ = -1;
    for (Vertex* v : this->inst_->vertices_) {

        this->bestrank_[v] = -1;
        this->bestnbrefs_[v] = -1;
    }
    this->objvalue_ = inst->nbvertices_;
    this->cliquecutsmaxsize_ = inst->nbvertices_;
}
//

DiscretizationSolver::DiscretizationSolver(Instance* inst, Algorithm algo, Problem pb, Model mod, std::string optionfile, float timelimit) :
problem_(pb), algo_(algo), modeltype_(mod), timelimit_(timelimit), inst_(inst) {
    this->bestnbfullref_ = -1;
    this->bestcliqueid_ = -1;
    for (Vertex* v : this->inst_->vertices_) {

        this->bestrank_[v] = -1;
        this->bestnbrefs_[v] = -1;
    }
    this->objvalue_ = inst->nbvertices_;
    if (!optionfile.empty())
    {
        std::ifstream infile(optionfile.c_str());
        if (!infile.is_open()) {
            throwError("revorder: error: the option file could not be opened");
        }
        // read the file line by line
        // one line contains the data relative to one option
        std::string line;
        while (std::getline(infile, line))
        {
            std::istringstream iss(line);
            char buf[100];
            std::string optionname;
            int optionvalue;
            
            //read the line and detect format errors
            if (!(iss >> buf)){
                std::cerr << "revorder: error: there is a mistake with the format of the option file" <<std::endl;
                throw;
            } // error
            optionname = std::string(buf);
            //
            if (optionname.find("cliquecutsmaxsize") !=std::string::npos) {
                if (!(iss >> optionvalue)){
                    std::cerr << "revorder: error: there is a mistake with the format of the option file" <<std::endl;
                    throw;
                } // error
                this->cliquecutsmaxsize_ = std::min(optionvalue, inst->nbvertices_);
                continue;
            }
            else if (optionname.find("initialcyclesize") !=std::string::npos) {
                if (!(iss >> optionvalue)){
                    std::cerr << "revorder: error: there is a mistake with the format of the option file" <<std::endl;
                    throw;
                } // error
                this->initialcyclesize_ = std::min(optionvalue, inst->nbvertices_);
                continue;
            }
            else if (optionname.find("decomposecomponents") !=std::string::npos) {
                if (!(iss >> optionvalue)){
                    std::cerr << "revorder: error: there is a mistake with the format of the option file" <<std::endl;
                    throw;
                } // error
                this->decomposecomponents_ = (bool) optionvalue;
                continue;
            }
        }
    }
    else {
        this->decomposecomponents_ = false;
        this->cliquecutsmaxsize_ = inst->nbvertices_;
        this->initialcyclesize_ = 3;
    }
}
//

DiscretizationSolver::~DiscretizationSolver() {

    for (Clique* clique : this->cliques_) delete clique;
    cliques_.clear();
}
//
// abstract solve method

int DiscretizationSolver::solve() {
    return 0;
}
//
// get the objective value of a given vertex order

int DiscretizationSolver::getobjvalue(std::map<Vertex*, int> rank, std::map<Vertex*, int> nbrefs) {

    int nbfullref = 0;
    for (Vertex* u : this->inst_->vertices_) {

        if (nbrefs[u] >= this->inst_->U()) nbfullref += 1;
    }
    return this->inst_->nbvertices_ - this->inst_->L() - nbfullref;

}
//
// get number of fully referenced vertices in an order

int DiscretizationSolver::getnbfullref(std::map<Vertex*, int> rank, std::map<Vertex*, int> nbrefs) {

    int nbfullref = 0;
    for (Vertex* u : this->inst_->vertices_) {

        if (nbrefs[u] >= this->inst_->U()) nbfullref += 1;
    }
    return nbfullref;

}
//
// return true if objval1 is better than objval2: this method is necessary, because the problem can be a minimization or a maximization

bool DiscretizationSolver::isabetterobjvalue(int objval1, int objval2) {

    return objval1 > objval2;
}
//
// compute the best greedy solution among those starting from the potential initial cliques and return true if the instance is feasible
// if allcliques is set to true, check every clique, otherwise, stop with the first feasible revorder

bool DiscretizationSolver::greedysolve(bool allcliques) {

    int U = this->inst_->U();
    int L = this->inst_->L();
    int nbvertices = this->inst_->nbvertices_;
    this->bestnbfullref_ = -1;
    this->bestcliqueid_ = -1;
    for (Vertex* u : this->inst_->vertices_) {
        this->bestisfullref_[u] = false;
        this->bestrank_[u] = -1;
        this->bestnbrefs_[u] = 0;
    }
    this->bestfullref_.clear();
    this->objvalue_ = this->inst_->nbvertices_;
    std::vector<Clique*> feasiblecliques;

    //
    // for each clique apply the constructive greedy that iteratively adds the
    // vertex with largest number of references
    std::map<Vertex*, int> rank;
    std::map<Vertex*, int> nbrefs;
    this->greedynbfullrefperclique_.clear();
    int feasiblecliqueid = -1;
    for (int initialcliqueindex = 0; initialcliqueindex < this->cliques_.size(); initialcliqueindex++) {
        //
        // data initialization
        Clique* initialclique = this->cliques_[initialcliqueindex];
        int currentrank = 1;
        std::vector<Vertex*> verticesleft = this->inst_->vertices_;
        int nbfullref = 0;
        for (Vertex * v : this->inst_->vertices_) {
            nbrefs[v] = 0;
            rank[v] = -1;
        }
        //
        // treat the initial clique first
        //
        // set the ranks of the vertices in the clique and set them as references of their neighbors
        for (Vertex* u : initialclique->vertices_) {
            rank[u] = currentrank++;
            for (Vertex* v : u->neighbors_) {
                if (rank[v] == -1) nbrefs[v]++;
            }
        }
        //
        // record the partially and fully referenced vertices
        std::vector<Vertex*> partialref;
        std::vector<Vertex*> fullref;
        for (Vertex* v : this->inst_->vertices_) {
            if (rank[v] == -1) {
                if (nbrefs[v] >= U) {
                    fullref.push_back(v);
                } else if (nbrefs[v] >= L) {
                    partialref.push_back(v);
                }
            }
        }
        //
        // include the vertex with largest number of references in the order until every vertex is treated or infeasibility of the initial clique is proved
        while (currentrank <= nbvertices) {
            if (!fullref.empty()) {
                Vertex* u = fullref.back();
                if (initialclique->initialpartialrefs_.empty()) {
                    initialclique->initialfullrefs_.push_back(u);
                }
                rank[u] = currentrank++;
                nbfullref++;
                fullref.pop_back();
                for (Vertex* v : u->neighbors_) {
                    if (rank[v] >= 1) continue;
                    nbrefs[v]++;
                    if (nbrefs[v] == U) {
                        fullref.push_back(v);
                    } else if (nbrefs[v] == L) {
                        partialref.push_back(v);
                    }
                }
            } else if (!partialref.empty()) {
                while (rank[partialref.back()] >= 1) {
                    partialref.pop_back();
                    if (partialref.empty()) {
                        break;
                    }
                }
                if (partialref.empty()) break;
                
                if (initialclique->initialpartialrefs_.empty()) {
                    for (Vertex* v: partialref) {
                        initialclique->initialpartialrefs_.push_back(v);
                    }
                }
                Vertex * u = partialref.back();
                rank[u] = currentrank++;
                partialref.pop_back();
                for (Vertex* v : u->neighbors_) {
                    if (rank[v] >= 0) continue;
                    nbrefs[v]++;
                    if (nbrefs[v] == U) {
                        fullref.push_back(v);
                    } else if (nbrefs[v] == L) {
                        partialref.push_back(v);
                    }
                }
            } else break;
        }

        if (currentrank - 1 == nbvertices) {
            feasiblecliqueid++;
            feasiblecliques.push_back(initialclique);
            this->greedynbfullrefperclique_.push_back(nbfullref);
            initialclique->greedyobjvalue_ = this->getobjvalue(rank, nbrefs);
            //
#ifdef VERBOSE
            std::cout << "revorder: greedy: clique computed solution with value " << initialclique->greedyobjvalue_ << std::endl;
#endif
            //
            // record the solution if a new incumbent is found
            if ((currentrank - 1 == nbvertices) && (initialclique->greedyobjvalue_ < this->objvalue_)) {
                //
#ifdef VERBOSE
                std::cout << "revorder: greedy: one new clique with cost " << initialclique->greedyobjvalue_ << '\n';
#endif
                //
                this->objvalue_ = initialclique->greedyobjvalue_;
                this->bestcliqueid_ = feasiblecliqueid;
                this->bestfullref_.clear();
                this->bestnbfullref_ = 0;
                for (Vertex* v : this->inst_->vertices_) {
                    this->bestrank_[v] = rank[v];
                    this->bestnbrefs_[v] = nbrefs[v];
                    if (this->bestnbrefs_[v] >= U) {
                        this->bestisfullref_[v] = true;
                        this->bestfullref_.push_back(v);
                        this->bestnbfullref_++;
                    } else {
                        this->bestisfullref_[v] = false;
                    }
                }
                //
                // if we do not need to check every clique, stop with the first feasible revorder
                if (!allcliques) break;
            }
        } else {
#ifdef VERBOSE
            std::cout << "revorder: greedy: this clique is not a feasible start" << std::endl;
#endif
        }
    }
    //
    // build the list of cliques that will be used in the optimization
    this->cliques_.clear();
    for (Clique* clique : feasiblecliques) {
        this->cliques_.push_back(clique);
    }
    //
    // update the list of cliques ids each vertex is a member of
    for (Vertex* v : this->inst_->vertices_) {
        this->cliqueslist_[v].clear();
        this->cliqueidslist_[v].clear();
    }
    int id = 0;
    for (Clique* clique : this->cliques_) {
        for (Vertex* v : clique->vertices_) {
            this->cliqueslist_[v].push_back(clique);
            this->cliqueidslist_[v].push_back(id);
        }
        id++;
    }
    std::cout << "revorder: greedy: number of cliques after greedy solution = " << this->cliques_.size() << std::endl;
    //
    // Stop here and return false if the problem is infeasible
    if (this->bestnbfullref_ == -1) {
        std::cout << "revorder: greedy: no solution was found with the greedy, the problem is infeasible" << std::endl;
        return false;
    }
    //
    // Record the solution found
    std::cout << "revorder: greedy: number of part. ref. vertices in the best solution = " << this->objvalue_ << std::endl;
    this->greedynbfullref_ = this->bestnbfullref_;

    return true;
}
//
//
// run the search for a discretization order with the input options

bool DiscretizationSolver::discretizationorder(float timelimit) {

    bool feas_status = true;

    if (this->algo_ == GREEDY) {
        //
        // enumerate the potential initial cliques
        this->inst_->computeneighbors();
        feas_status = this->enumeratecliques(this->inst_->L() + 1);
        if (feas_status == false) {
            this->isfeasible_ = false;
            return false;
        }
        this->eliminateredundantcliques();
        //
        // run the greedy search
        feas_status = greedysolve();
        this->isfeasible_ = feas_status;
        return feas_status;
    }


    switch (this->algo_) {
        case IP:
            feas_status = minpartialip(timelimit);
            break;
        case CP:
            feas_status = cpvertex(timelimit);
            break;

        default:
            std::cout << "revorder: error with the problem and/or algorithm fed to the solver" << std::endl;
            throw;
    }


    return feas_status;
}
//
// greedy algorithm to get a reference primal bound
// try all the possible initial cliques as starting points of the greedy

bool DiscretizationSolver::greedymipstart(IloCplex & cplex) {
    //
    // set the best greedy solution as primal solution for the mip solution
    IloNumArray mipstartvalues(cplex.getEnv());
    IloNumVarArray mipstartvariables(cplex.getEnv());
    if (this->modeltype_ == VERTEXRANK) {
        for (Vertex* v : this->inst_->vertices_) {
            for (int k = 0; k < this->inst_->nbvertices_; k++) {
                mipstartvariables.add(this->hasrank_[v][k]);
                if (this->bestrank_[v] == k + 1) mipstartvalues.add(1);
                else mipstartvalues.add(0);
            }
        }
    } else {
        if (this->modeltype_ == WITNESS) {
            std::map<Vertex*, int> nbwitness;
            for (Vertex* u : this->inst_->vertices_) nbwitness[u] = 0;

            for (Vertex* u : this->inst_->vertices_) {
                if (this->bestrank_[u] <= this->inst_->L() + 1) {
                    for (Vertex* v : u->neighbors_) {
                        mipstartvariables.add(this->isbefore_[u][v]);
                        mipstartvalues.add(1);
                        nbwitness[v] += 1;
                    }
                }
            }
            for (Vertex* u : this->inst_->vertices_) {
                if (this->bestrank_[u] >= this->inst_->L() + 2) {
                    for (Vertex* v : u->neighbors_) {
                        mipstartvariables.add(this->isbefore_[u][v]);
                        if ((this->bestrank_[u] < this->bestrank_[v])) {
                            if ((this->bestisfullref_[v]) && (nbwitness[v] < this->inst_->U())) {
                                mipstartvalues.add(1);
                                nbwitness[v] += 1;
                            } else if ((!this->bestisfullref_[v]) && (nbwitness[v] < this->inst_->L())) {
                                mipstartvalues.add(1);
                                nbwitness[v] += 1;
                            } else mipstartvalues.add(0);
                        } else mipstartvalues.add(0);
                    }
                }
            }
        } else {
            for (Vertex* u : this->inst_->vertices_) {
                for (Vertex* v : u->neighbors_) {
                    if (u->id_ >= v->id_) continue;
                    mipstartvariables.add(this->isbefore_[u][v]);
                    mipstartvariables.add(this->isbefore_[v][u]);
                    if (this->bestrank_[u] < this->bestrank_[v]) {
                        mipstartvalues.add(1);
                        mipstartvalues.add(0);
                    } else {
                        mipstartvalues.add(0);
                        mipstartvalues.add(1);
                    }
                }
            }
        }
        //
        // set the clique variables and corresponding clique variables
        if (this->modeltype_ == WITNESS) {
            for (Vertex* v : this->inst_->vertices_) {
                mipstartvariables.add(this->isvertexinclique_[v]);
                if (this->bestrank_[v] <= this->inst_->L() + 1) mipstartvalues.add(1);
                else mipstartvalues.add(0);
            }
        } else {
            for (int c = 0; c < this->cliques_.size(); c++) {
                mipstartvariables.add(this->isclique_[c]);
                if (c == this->bestcliqueid_) mipstartvalues.add(1);
                else mipstartvalues.add(0);
            }
        }
        //
        // full-ref variables
        for (Vertex* u : this->inst_->vertices_) {
            mipstartvariables.add(this->isfullref_[u]);
            if ((this->bestisfullref_[u]) || (this->bestrank_[u] <= this->inst_->L() + 1)) mipstartvalues.add(1);
            else mipstartvalues.add(0);
        }
        //
        // rank variables for RANKS model
        for (Vertex* v : this->inst_->vertices_) {
            if (this->modeltype_ == RANKS) {
                for (int k = 0; k < this->inst_->nbvertices_; k++) {
                    mipstartvariables.add(this->hasrank_[v][k]);
                    if (this->bestrank_[v] == k + 1) mipstartvalues.add(1);
                    else mipstartvalues.add(0);
                }
                mipstartvariables.add(this->rank_[v]);
                mipstartvalues.add(this->bestrank_[v] - 1);
            }
        }
    }

    cplex.addMIPStart(mipstartvariables, mipstartvalues);
    mipstartvariables.end();
    mipstartvalues.end();

    return true;
}


//
// contraint programming CP^VERTEX from Bodur and MacNeil, 2019

bool DiscretizationSolver::cpvertex(float timelimit) {

    char name[50];
    int nbvertices = this->inst_->nbvertices_;
    int L = this->inst_->L();
    int U = this->inst_->U();
    IloEnv env;
    IloTimer cpuClockTotal(env);
    IloTimer cpuClockInit(env);
    cpuClockTotal.start();
    cpuClockInit.start();


    std::cout << std::endl;
    std::cout << "revorder: ----------------------------------------------------------" << std::endl;
    std::cout << "\nrevorder: cp: Run preprocessing procedures " << std::endl << std::endl;

    //
    // compute adjacency lists
    this->inst_->computeneighbors();
    //
    // identify the set of feasible initial cliques
    if (!this->enumeratecliques(this->inst_->L() + 1)) {
        this->isfeasible_ = false;
        return false;
    }
    this->eliminateredundantcliques();
    //
    // run the greedy algorithm for initial solution and reduction of initial cliques
    if (!this->greedysolve()) {
        std::cout << "Greedy solve found the problem to be infeasible at the time of mip warm start." << std::endl;
        return false;
    }

    //
    // initialize the model
    IloModel model(env);
    //
    // assign one unique index in {0,...,n-1} to each vertex
    std::map<Vertex*, int> indv;
    int ind = 0;
    for (Vertex* u : this->inst_->vertices_) {
        indv[u] = ind++;
    }

    //
    // initialize adjacency matrix
    IloIntArray adjmatrix(env, nbvertices * nbvertices);
    for (Vertex* u : this->inst_->vertices_) {
        for (Vertex* v : this->inst_->vertices_) {
            if (indv[v] < indv[u]) continue;
            if (u->isneighbor_[v]) {
                adjmatrix[nbvertices * indv[u] + indv[v]] = 1;
                adjmatrix[nbvertices * indv[v] + indv[u]] = 1;
            } else {
                adjmatrix[nbvertices * indv[u] + indv[v]] = 0;
                adjmatrix[nbvertices * indv[v] + indv[u]] = 0;
            }
        }
    }
    //
    // variables that set the rank of each vertex in the order
    //
    // isrankfullref[k] =1 if vertex v is fully-referenced and 0 otherwise
    IloBoolVarArray isrankfullref(env, nbvertices);
    for (int k = 0; k < nbvertices; k++) {
        isrankfullref[k] = IloBoolVar(env);
        sprintf(name, "isrankfullref_%i", k);
        isrankfullref[k].setName(name);
    }
    //
    // indatrank[r] = index of the vertex at rank r
    IloIntVarArray indatrank(env, nbvertices);
    for (int k = 0; k < nbvertices; k++) {
        indatrank[k] = IloIntVar(env, 0, nbvertices - 1);
        sprintf(name, "indatrank_%i", k);
        indatrank[k].setName(name);
    }
    //
    // objective function
    IloExpr obj(env);
    obj += 1;
    for (int k = L + 1; k < nbvertices; k++) {
        obj += (1 - isrankfullref[k]);
    }
    model.add(IloMinimize(env, obj));

    for (int i = 0; i <= L - 1; i++) {
        for (int j = i + 1; j <= L; j++) {
            model.add(IloElement(adjmatrix, nbvertices * indatrank[i] + indatrank[j]) == 1);
        }
    }

    for (int r = L + 1; r < nbvertices; r++) {
        IloExpr sumrefs(env);
        for (int i = 0; i <= r - 1; i++) {
            sumrefs += IloElement(adjmatrix, nbvertices * indatrank[i] + indatrank[r]);
        }
        model.add(sumrefs - (U - L) * isrankfullref[r] >= L);
        sumrefs.end();
    }

    model.add(IloAllDiff(env, indatrank));
    //
    // set as not in clique, all the vertices that do not belong to one of the possible clique
    std::map < Vertex*, bool> okforclique;
    for (Vertex* v : this->inst_->vertices_) okforclique[v] = false;
    for (Clique* c : this->cliques_) {
        for (Vertex* v : c->vertices_) {
            okforclique[v] = true;
        }
    }
    for (Vertex* v : this->inst_->vertices_) {
        if (!okforclique[v]) {
            for (int r = 0; r <= L; r++) {
                model.add(indatrank[r] != indv[v]);
            }
        }
    }
    //
    // load model in the CP solver
    IloCP cp(model);
    //
    // provide starting point with greedy solution
    //
    IloSolution sol(env);
    for (Vertex* v : this->inst_->vertices_) {
        if (this->bestrank_[v] - 1 >= L + 1) {
            sol.setValue(isrankfullref[this->bestrank_[v] - 1], bestisfullref_[v]);
        }
        sol.setValue(indatrank[this->bestrank_[v] - 1], indv[v]);
    }
    cp.setStartingPoint(sol);

    cpuClockInit.stop();
    // 
    // set CP parameter
    cp.setParameter(IloCP::LogVerbosity, IloCP::Normal);
    cp.setParameter(IloCP::Workers, 1);
    cp.setParameter(IloCP::TimeLimit, timelimit - cpuClockInit.getTime());
#ifdef VERBOSE
    cp.setParameter(IloCP::LogVerbosity, IloCP::Normal);
#endif
    IloTimer cpuClock(env);
    cpuClock.start();
    try {
        cp.solve();
    } catch (IloException& e) {
        std::cerr << e << std::endl;
    }
    cpuClockTotal.stop();

    IloAlgorithm::Status status = cp.getStatus();
    std::cout << "revorder: ----------------------------------------------------------\n" << std::endl;
    std::cout << "revorder: CP optimizer solution status = " << status << std::endl;
    isfeasible_ = (status == IloAlgorithm::Feasible) || (status == IloAlgorithm::Optimal);
    isoptimal_ = (status == IloAlgorithm::Optimal);
    if (isfeasible_) {
        //
        // retrieve the information on the solution process
        this->objvalue_ = cp.getObjValue();
        totaltime_ = cpuClock.getTime();
        treatednodes_ = cp.getInfo(IloCP::NumberOfBranches);

        //
        // get the optimal order
        for (int r = 0; r < this->inst_->nbvertices_; r++) {
            this->bestrank_[this->inst_->vertices_[cp.getValue(indatrank[r])]] = r;
        }
        //
        // get the number of references of each vertex
        for (Vertex* u : this->inst_->vertices_) {
            this->bestnbrefs_[u] = 0;
            for (Vertex* v : u->neighbors_) {
                if (this->bestrank_[v] < this->bestrank_[u]) this->bestnbrefs_[u]++;
            }
        }
        //
        // get fully referenced vertices
        this->bestfullref_.clear();
        this->bestnbfullref_ = 0;
        this->bestisfullref_.clear();
        for (Vertex* u : this->inst_->vertices_) {
            if (this->bestnbrefs_[u] >= this->inst_->U()) {
                this->bestisfullref_[u] = true;
                this->bestfullref_.push_back(u);
                this->bestnbfullref_++;
            } else this->bestisfullref_[u] = false;
        }
        //
        // summary of cp execution
        std::cout << "revorder: total cpu time                    = " << cpuClockTotal.getTime() << " s" << std::endl;
        std::cout << "revorder: initialization cpu time           = " << cpuClockInit.getTime() << " s" << std::endl;
        std::cout << "revorder: number of branches explored       = " << treatednodes_ << std::endl;
        std::cout << "revorder: value of the objective function   = " << objvalue_ << std::endl;
        std::cout << "revorder: verification: number of part. ref. vertices = " << this->inst_->nbvertices_ - this->inst_->L() - this->bestnbfullref_ << std::endl;
    } else {
        std::string fileInfeasible("InfeasibleModel.lp");
        cp.exportModel(fileInfeasible.c_str());
        std::cout << "revorder: no feasible solution was found during the optimization!" << std::endl;
        std::cout << "revorder: check the file " << fileInfeasible << " to see the corresponding model" << std::endl;
    }
    //
    // delete CP objects
    cp.end();
    model.end();
    env.end();

    return isfeasible_;
}
//
// create the cplex model for the IP formulation described in Bodur an MacNeil, 2019

void DiscretizationSolver::defineminpartial_vertexrank(IloModel & model) {
    IloEnv env = model.getEnv();
    char name[50];
    int nbvertices = this->inst_->nbvertices_;
    int L = this->inst_->L();
    int U = this->inst_->U();

    //////////////////////////////////////////////////////////////////////////////
    // Declare the variables
    //////////////////////////////////////////////////////////////////////////////
    //
    // variables that set the rank of each vertex in the order
    // - hasrank_[v,k] = 1 if vertex i has rank k in the discretization order
    //  in this model, we use these variables only for the first L vertices
    //  of the cliques (if modeltype <= 1)
    for (Vertex* v : this->inst_->vertices_) {
        hasrank_[v] = IloBoolVarArray(env, nbvertices);
        for (int k = 0; k < nbvertices; k++) {
            sprintf(name, "hasrank_%i_%i", v->id_, k);
            hasrank_[v][k].setName(name);
        }
    }
    //
    // isrankfullref[k] =1 if vertex v is fully-referenced and 0 otherwise
    IloBoolVarArray isrankfullref(env, nbvertices);
    for (int k = 0; k < nbvertices; k++) {
        isrankfullref[k] = IloBoolVar(env);
        sprintf(name, "isrankfullref_%i", k);
        isrankfullref[k].setName(name);
    }
    //
    // isfullrefatrank[u,k] = 1 if vertex u is fully-referenced and is at rank k
    BoolVarMapArray isfullrefatrank;
    for (Vertex* v : this->inst_->vertices_) {
        isfullrefatrank[v] = IloBoolVarArray(env, nbvertices);
        for (int k = 0; k < this->inst_->nbvertices_; k++) {
            sprintf(name, "isfullrefatrank%i_%i", v->id_, k);
            isfullrefatrank[v][k].setName(name);
        }
    }

    //////////////////////////////////////////////////////////////////////////////
    // Set the cplex model that will be solved
    //////////////////////////////////////////////////////////////////////////////
    //
    // objective function
    IloExpr obj(env);
    obj += 1;
    for (int k = 0; k < nbvertices; k++) {
        obj += (1 - isrankfullref[k]);
    }
    model.add(IloMinimize(env, obj));

    //
    // each rank of the order is taken by exactly one vertex
    IloRangeArray ctaryOneVertexPerRank(env);
    for (int k = 0; k < nbvertices; k++) {
        IloExpr sumvertices(env);
        for (Vertex* v : this->inst_->vertices_) sumvertices += hasrank_[v][k];
        ctaryOneVertexPerRank.add(sumvertices == 1);
        sprintf(name, "ctaryOneVertexPerRank_%i", k);
        ctaryOneVertexPerRank[k].setName(name);
    }
    model.add(ctaryOneVertexPerRank);
    //
    // each vertex has exactly one rank and his rank is computed
    IloRangeArray ctaryOneRankPerVertex(env);
    IloRangeArray ctaryComputeRank(env);
    int nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
        IloExpr sumhasrank(env);
        for (int k = 0; k < nbvertices; k++) sumhasrank += hasrank_[u][k];
        ctaryOneRankPerVertex.add(sumhasrank == 1);
        sprintf(name, "ctaryOneRankPerVertex_%i", u->id_);
        ctaryOneRankPerVertex[nbcons].setName(name);
    }
    model.add(ctaryOneRankPerVertex);
    //
    // Guarantee that isfullrefatrank[v][k] =1 only if v has rank k and has at least U references
    IloRangeArray ctaryIsFullRefAtLevel(env);
    for (int k = 0; k <= L; k++) {
        ctaryIsFullRefAtLevel.add(isrankfullref[k] == 1);
    }
    for (int k = 0; k < this->inst_->nbvertices_; k++) {
        for (Vertex* v : this->inst_->vertices_) {
            ctaryIsFullRefAtLevel.add(hasrank_[v][k] - (1 - isrankfullref[k]) - isfullrefatrank[v][k] <= 0);
        }
    }
    model.add(ctaryIsFullRefAtLevel);
    //
    // Discretization constraints: every vertex must be partially-referenced and those with more than U references are fully-referenced
    IloRangeArray ctaryLrefs(env);
    IloRangeArray ctaryClique(env);
    IloRangeArray ctaryFullRef(env);
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
        //
        // first deal with the vertices of the initial clique
        for (int k = 0; k <= L; k++) {
            IloExpr sumrefs(env);
            for (Vertex* v : u->neighbors_) {
                for (int l = 0; l <= k - 1; l++) sumrefs += hasrank_[v][l];
            }
            ctaryClique.add(sumrefs - k * hasrank_[u][k] >= 0);
        }
        //
        // then deal with the rest of the order
        for (int k = L + 1; k < nbvertices; k++) {
            IloExpr sumrefs(env);
            for (Vertex* v : u->neighbors_) {
                for (int l = 0; l <= k - 1; l++) sumrefs += hasrank_[v][l];
            }
            ctaryLrefs.add(sumrefs - L * hasrank_[u][k] >= 0);
            ctaryFullRef.add(sumrefs - U * isfullrefatrank[u][k] >= 0);
        }
    }
    model.add(ctaryClique);
    model.add(ctaryLrefs);
    model.add(ctaryFullRef);
    //
    // set as not in clique, all the vertices that do not belong to one of the possible clique