MILP_Model.h

#ifndef MILP_MODEL_H
#define MILP_MODEL_H

#include <algorithm>
#include <cassert>
#include <fstream>
#include <map>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
#include <string.h>
#include <unistd.h>
#include <regex>

using namespace std;

class Milp_Model
{
public:

    enum VarType {REAL, INTEGER, BOOLEAN};  
    enum RowType {LEQ, GEQ, EQ};            
    enum Status {OPTIMAL, NONOPTIMAL, UNFEASIBLE, UNBOUNDED, UNKNOWN, ERROR};   

    //SHAB_q: I don't understand what 'Term' means.
    //SHAB_ans: the pair literally mean: coefficient * vecVars[var index]. like 2a, 3x, etc.
    using Term = pair<double, int>;         
    using vecTerms = vector<Term>;          
    using vecVars = vector<int>;            

    Milp_Model(const string& solver="") {
        init(solver);
    }

    const string& getError() const {
        return errorMsg;
    }

    void setError(const string& err = "") {
        errorMsg = err;
    }

    void setMinimize() {
        MinMax = true;
    }

    void setMaximize() {
        MinMax = false;
    }

    void setEpsilon(double eps) {
        epsilon = eps;
    }

    int numVariables() const {
        return Vars.size();
    }

    int numConstraints() const {
        return Matrix.size();
    }

    int newRealVar(string name="",
                   double lower_bound = 0.0,
                   double upper_bound = -1.0) {
        return newVar(name, REAL, lower_bound, upper_bound);
    }

    int newIntegerVar(string name="",
                      double lower_bound = 0,
                      double upper_bound = -1) {
        return newVar(name, INTEGER, lower_bound, upper_bound);
    }

    int newBooleanVar(string name="") {
        return newVar(name, BOOLEAN, 0, 1);
    }

    Status getStatus() const {
        return stat;
    }

    double getObj() const {
        return obj;
    }

    double operator[](int i) const {
        assert (i >= 0 and i < Vars.size());
        return Vars[i].value;
    }

    double operator[](const string& name) const {
        assert (Name2Var.count(name) > 0);
        return Vars[Name2Var.at(name)].value;
    }

    bool isTrue(int i) const {
        assert (i >= 0 and i < Vars.size() and Vars[i].type == BOOLEAN);
        return Vars[i].value > 0.5;
    }

    bool isFalse(int i) const {
        return not isTrue(i);
    }

    int newRow(char type = '<', double rhs = 0.0, const string& name = "") {
        RowType t;
        if (type == '<') t = LEQ;
        else if (type == '>') t = GEQ;
        else if (type == '=') t = EQ;
        else assert(false);
        Matrix.push_back(Row {name, t, rhs});
        return Matrix.size() - 1;
    }

    int newRow(const vecTerms& terms, char type = '<', double rhs = 0.0, const string& name = "") {
        int r = newRow(type, rhs, name);
        for (auto t: terms) newTerm(r, t.first, t.second);
        return r;
    }

    void newTerm(int rowIndex, double coeff, int varIndex) {
        assert (rowIndex >= 0);
        assert (rowIndex < Matrix.size());
        assert (varIndex >= 0);
        assert (varIndex < Vars.size());
        assert (rowIndex >= 0 and rowIndex < Matrix.size() and
                varIndex >= 0 and varIndex < Vars.size());
        Matrix[rowIndex].vecRow.push_back( {coeff, varIndex});
    }

    void setRHS(int rowIndex, double rhs) {
        assert (rowIndex >= 0 and rowIndex < Matrix.size());
        Matrix[rowIndex].rhs = rhs;
    }

    int atMostOne(const vecVars& vars, const string& name="") {
        int r = newRow('<', 1.0, name);
        for (int v: vars) newTerm(r, 1.0, v);
        return r;
    }

    int atLeastOne(const vecVars& vars, const string& name="") {
        int r = newRow('>', 1.0, name);
        for (int v: vars) newTerm(r, 1.0, v);
        return r;
    }

    int exactlyOne(const vecVars& vars, const string& name="") {
        int r = newRow('=', 1.0, name);
        for (int v: vars) newTerm(r, 1.0, v);
        return r;
    }

    int equalSum(const vecVars& x, const vecVars& y, const string& name="") {
        int r = newRow('=', 0.0, name);
        for (int v: x) newTerm(r, 1.0, v);
        for (int v: y) newTerm(r, -1.0, v);
        return r;
    }

    int allEqual(const vecVars& v) {
        assert (v.size() > 1);
        int row = newRow( {{1, v[0]}, {-1, v[1]}}, '=');
        for (unsigned int i = 2; i < v.size(); ++i) newRow( {{1, v[0]}, {-1, v[i]}}, '=');
        return row;
    }

    int implies(int x, const vecVars& y) {
        int r = newRow('>', 0.0);
        newTerm(r, -1.0, x);
        for (int v: y) newTerm(r, 1.0, v);
        return r;
    }

    int implies(const vecVars& x, int y) {
        int first = -1;
        for (int v: x) {
            int r = newRow('>', 0.0);
            if (first < 0) first = r;
            newTerm(r, -1.0, v);
            newTerm(r, 1.0, y);
        }
        return first;
    }


    void newCostTerm(double coeff, int varIndex) {
        assert (varIndex >= 0 and varIndex < Vars.size());
        Cost.push_back( {coeff, varIndex});
    }

    bool writeLP(const string& filename) {
        ofstream f;
        f.open(filename);
        if (not f.is_open()) {
            setError("Could not open file " + filename + ".");
            return false;
        }
        writeLP(f);
        f.close();
        return true;
    }

    bool solve(int timelimit = -1) {

        if (solver != "cbc" and solver != "glpsol" and solver != "gurobi_cl") { //Carmine 25.02.22 gurobi_cl included as MILP solver
            setError("Unkonwn solver " + solver + ".");
            return false;
        }

        string exec = "which " + solver + " >/dev/null 2>&1";
        if (system(exec.c_str())) {
            setError("Solver " + solver + " is not available.");
            return false;
        }

        // The solver seems to be working. Let us solve the MILP model.
        string lpfile = createTempFilename("MILP_lpmodel", ".lp");
        string outfile;
        if(solver == "gurobi_cl")  //Carmine 25.02.22 gurobi_cl accepts solution in sol format
            outfile = createTempFilename("MILP_solution", ".sol");
        else
            outfile = createTempFilename("MILP_solution", ".gsol");
        string command = writeCommand(solver, lpfile, outfile, timelimit) + " >/dev/null 2>&1";
        writeLP(lpfile);

        int status = system(command.c_str());

        if (status != 0) {
            deleteTempFilename(outfile);
            setError("Error when executing " + solver + ".");
            return false;
        }
        bool opt_status; //Carmine 25.02.22 Important to check if the status optimization is OPTIMAL or SUBOPTIMAL
        if (solver == "cbc") opt_status = readCbcSolution(outfile);
        else if (solver == "glpsol") opt_status = readGlpsolSolution(outfile);
        else if (solver == "gurobi_cl") opt_status = readGurobiSolution(outfile); //Carmine 25.02.22 gurobi_cl included as MILP solver
        else assert(false);

        if(!opt_status)
            cout << "*ERROR* MILP solution is UNFEASIBLE or UNBOUNDED" << endl;

        //deleteTempFilename(lpfile);
        //deleteTempFilename(outfile);
        return true;
    }

    bool init(const string& solver = "") {
        MinMax=true;
        epsilon=10e-10;
        numRealVars=numIntegerVars=numBooleanVars=0;
        numEmptyRows = 0;
        numUsedVars = 0;
        stat=UNKNOWN;
        errorMsg = "";
        Cost.clear();
        Vars.clear();
        Matrix.clear();
        Name2Var.clear();
        Name2delays.clear();
        return find_solver(solver);
    }

    void writeOutDelays(const string& file_name) {
        ofstream file;
        file.open(file_name);
        for (auto const& values: Name2delays){
            file << values.second << "  " <<  Vars[values.first].value << endl;
        }
        file.close();

        file.open("tmp_delays.txt");
        int i, ord = Vars.size();
        for (i = 0; i< ord; i++){
            file <<  Vars[i].name << "  " <<  Vars[i].value << endl;
        }
        file.close();
    }

    string getVarName(int id) {
        return Vars[id].name;
    }

private:

    struct Var {
        string name;
        VarType type;
        double lower_bound;
        double upper_bound; // unbounded if upper_bound < lower_bound
        double value;
    };

    struct Row {
        string name;
        RowType type;           
        double rhs;             
        vecTerms vecRow;    
    };

    string solver;      
    bool MinMax;        
    vecTerms Cost;      
    vector<Var> Vars;   
    vector<Row> Matrix; 
    int numEmptyRows;   
    vector<int>appearanceOrder; 
    vector<bool>appeared;       
    int numUsedVars;            
    double epsilon;     
    int numRealVars;
    int numIntegerVars;
    int numBooleanVars;
    map<string, int> Name2Var;  
    Status stat;        
    double obj;        
    string errorMsg;    

    map<int, string> Name2delays; //Carmine 07.02.2022 map containing the output delay of the blocks

    bool find_solver(const string& s) {

        if (not s.empty() and s != "cbc" and s != "glpsol" and s != "gurobi_cl") { //Carmine 25.02.22 gurobi_cl included as MILP solver
            setError("Unkonwn solver " + solver + ".");
            return false;
        }
        
        // First try cbc
        if (s.empty() or s == "cbc") {
            string exec = "which cbc >/dev/null 2>&1";
            if (system(exec.c_str()) == 0) {
                solver = "cbc";
                return true;
            }

            if (s == "cbc") {
                setError("MILP solver cbc not found.");
                return false;
            }
        }

        // Next try glpsol
        if (s.empty() or s == "glpsol") {
            string exec = "which glpsol >/dev/null 2>&1";
            if (system(exec.c_str()) == 0) {
                solver = "glpsol";
                return true;
            }

            if (s == "glpsol") {
                setError("MILP solver glpsol not found.");
                return false;
            }
        }

        // Next try gurobi_cl
        if (s.empty() or s == "gurobi_cl") { //Carmine 25.02.22 gurobi_cl included as MILP solver
            string exec = "which gurobi_cl >/dev/null 2>&1";
            if (system(exec.c_str()) == 0) {
                solver = "gurobi_cl";
                return true;
            }

            if (s == "gurobi_cl") {
                setError("MILP solver gurobi_cl not found.");
                return false;
            }
        }

        // No more solvers for the moment
        setError("No MILP solver found.");
        return false;
    }

    int newVar(const string& name, VarType type, double lower_bound, double upper_bound) {
        string n = name.empty() ? "x" + to_string(Vars.size()) : name;
        if (Name2Var.count(n) > 0) {
            setError("Variable " + n + " multiply defined.");
            return -1;
        }
        Name2Var[n] = Vars.size();
        Vars.push_back(Var {n, type, lower_bound, upper_bound, 0});
        if (type == REAL) numRealVars++;
        else if (type == INTEGER) numIntegerVars++;
        else numBooleanVars++;
        return Vars.size() - 1;
    }

    int newOutDelay(int ind, string output) {
        if (Name2delays.count(ind) > 0)
            return 0;
        Name2delays[ind] = output;
        return 1;
    }

    void normalize() {
        normalizeRow(Cost);
        for (Row& r: Matrix) normalizeRow(r.vecRow);
    }

    void normalizeRow(vecTerms& r) {
        // Sorts the terms by variable index
        sort(r.begin(), r.end(), [](Term& t1, Term& t2) {
            return t1.second < t2.second;
        });

        // Now check consecutive terms and accumulate coefficients
        int i = 0;
        int last = r.size();
        while (i < last - 1) {
            int var = r[i].second;
            unsigned int k = i + 1;
            while (k < last and r[k].second == var) {
                r[i].first += r[k].first;
                r[k].first = 0;
                ++k;
            }
            i = k;
        }

        // Now remove the terms with zero coefficient
        unsigned int k = 0;  // to point at the first available slot
        for (i = 0; i < last; ++i) {
            if (abs(r[i].first) >= epsilon) r[k++] = r[i];
        }
        r.resize(k);
    }

    void writeTerms(ofstream& f, const vecTerms& terms) {
        assert (terms.size() > 0);
        double coeff = terms[0].first;
        if (abs(coeff) != 1.0) f << coeff << ' ';
        else if (coeff < 0) f << '-';
        int idx = terms[0].second;
        f << Vars[idx].name;
        if (not appeared[idx]) {
            appearanceOrder.push_back(idx);
            appeared[idx] = true;
            regex regexp ("timePath_(.)+_out[0-9]+");  //Carmine 07.02.2022 Extracting timing output pins indexes
            if (regex_search(Vars[idx].name, regexp)){
                newOutDelay(idx, Vars[idx].name);
            }
        }

        for (unsigned int i = 1; i < terms.size(); ++i) {
            coeff = terms[i].first;
            f << ' ' << (coeff < 0 ? '-' : '+') << ' ';
            if (abs(coeff) != 1) f << abs(coeff) << ' ';
            idx = terms[i].second;
            f << Vars[idx].name;
            if (not appeared[idx]) {
                appearanceOrder.push_back(idx);
                appeared[idx] = true;
                regex regexp ("timePath_(.)+_out[0-9]+");  //Carmine 07.02.2022 Extracting timing output pins indexes
                if (regex_search(Vars[idx].name, regexp)){
                    newOutDelay(idx, Vars[idx].name);
            }
            }
        }
    }

    void writeRow(ofstream& f, const Row& r) {
        if (r.vecRow.empty()) {
            ++numEmptyRows;
            return;
        }
        f << "  ";
        if (not r.name.empty()) f << r.name << ": ";
        writeTerms(f, r.vecRow);
        if (r.type == EQ) f << " = ";
        else if (r.type == GEQ) f << " >= ";
        else f << " <= ";
        f << r.rhs << endl;
    }

    void writeLP(ofstream& f) {
        normalize();
        appearanceOrder.clear();
        appeared = vector<bool>(Vars.size(), false);

        // Cost function: if no cost, create a fake cost function
        if (Cost.empty()) newCostTerm(0, 0);

        // Cost function
        f << (MinMax ? "Minimize" : "Maximize") << endl;
        f << "  ";
        writeTerms(f, Cost);
        f << endl;

        f << "Subject to" << endl;
        for (const Row& r: Matrix) writeRow(f, r);

        numUsedVars = 0;
        for (bool b: appeared) if (b) ++numUsedVars;

        bool need_bounds = false;
        for (const Var& v: Vars) {
            if (v.type != BOOLEAN and v.lower_bound <= v.upper_bound) {
                need_bounds = true;
                break;
            }
        }

        if (need_bounds) {
            f << "Bounds" << endl;
            for (const Var& v: Vars) {
                if (v.type != BOOLEAN and v.lower_bound <= v.upper_bound) {
                    f << "  " << v.lower_bound << " <= " << v.name << " <= " << v.upper_bound << endl;
                }
            }
        }

        if (numIntegerVars > 0) {
            f << "General" << endl << ' ';
            for (const Var& v: Vars) {
                if (v.type == INTEGER) f << ' ' << v.name;
            }
            f << endl;
        }

        if (numBooleanVars > 0) {
            f << "Binary" << endl << ' ';
            for (const Var& v: Vars) {
                if (v.type == BOOLEAN) f << ' ' << v.name;
            }
            f << endl;
        }

        f << "End" << endl;
    }

    static int readLine(ifstream& f, vector<string>& items) {
        string line;
        getline(f, line);
        istringstream myStream(line);
        items.clear();
        string s;
        while (myStream >> s) items.push_back(s);
        return items.size();
    }

    bool readCbcSolution(const string& filename) {
        ifstream f;
        f.open(filename);
        if (not f.is_open()) return false;

        vector<string> values;
        int nv = readLine(f, values);
        assert(nv == 2);

        int nrows = stoi(values[0]);
        int ncols = stoi(values[1]);

        assert (nrows == Matrix.size() - numEmptyRows + 1);
        assert (ncols == numUsedVars);

        // This is a tricky function that has to deal with the
        // fact that cbc sometimes puts three values per line
        // and sometimes only one (possibly when MILP model).
        nv = readLine(f, values);

        int st = stoi(values[0]);
        switch (st) {
        case 1:
            stat = UNBOUNDED;
            break;
        case 2:
            stat = NONOPTIMAL;
            break;
        case 4:
            stat = UNFEASIBLE;
            break;
        case 5:
            stat = OPTIMAL;
            break;
        default:
            assert(false);
        }

        if (stat != OPTIMAL and stat != NONOPTIMAL) {
            f.close();
            return false; //Carmine 25.02.22 if it is unfeasible or unbounded it should stop dynamatic execution
        }

        obj = stod(values.back());
        if (not MinMax) obj = -obj;

        // Skip the value of the rows
        for (unsigned int i = 0; i < nrows; ++i) readLine(f, values);

        // Read the variables
        for (unsigned int i = 0; i < ncols; ++i) {
            nv = readLine(f, values);
            int loc = nv > 1 ? 1 : 0;
            Vars[appearanceOrder[i]].value = stod(values[loc]);
        }

        f.close();
        return true;
    }

    bool readGlpsolSolution(const string& filename) {
        ifstream f;
        f.open(filename);
        if (not f.is_open()) return false;

        bool cost_read = false;
        vector<string> line;

        while (readLine(f, line) > 0) {
            char c = line[0][0];
            if (c != 's' and c != 'j') {
                // Skip line
                continue;
            } else if (c == 'j') {
                // Now we can read the value of a variable
                int nvar = stoi(line[1]);
                int loc = line.size() == 3 ? 2 : 3;
                Vars[appearanceOrder[nvar-1]].value = stod(line[loc]);
            } else { // c == 's'
                assert (not cost_read);
                cost_read = true;
                // Read the cost function
                int nrows = stoi(line[2]);
                int ncols = stoi(line[3]);
                char st = line[4][0];
                obj = stod(line.back());

                assert (nrows == Matrix.size() - numEmptyRows);
                assert (ncols == numUsedVars);

                switch (st) {
                case 'u':
                    stat = UNBOUNDED;
                    break;
                case 'f':
                    stat = NONOPTIMAL;
                    break;
                case 'n':
                    stat = UNFEASIBLE;
                    break;
                case 'o':
                    stat = OPTIMAL;
                    break;
                default:
                    assert(false);
                }
            }
        }

        f.close();
        return true;
    }

    bool readGurobiSolution(const string& filename) {
        ifstream f;
        f.open(filename);
        if (not f.is_open()) return false;

        vector<string> values;
        readLine(f, values); // first line only tells the optimization function cost of the solution

        string line;
        int numVars = 0;
        while(getline(f, line)){
            istringstream myStream(line);
            values.clear();
            string s;
            while (myStream >> s) values.push_back(s);
            Vars[Name2Var[values[0]]].value = stod(values[1]);
            numVars++;
        }
        
        //assert (numVars == numUsedVars);

        //TODO: check the optimization status of gurobi
        stat = OPTIMAL;
        
        f.close();
        return true;
    }

    string writeCommand(const string& solver, const string& lpfile, const string& solfile,
                        int timeout = -1) {
        ostringstream command;
        command << solver << ' ';
        if (solver == "cbc") {
            command << lpfile;
            if (timeout > 0) command << " sec " << timeout;
            command << " solve gsolution " << solfile;
        } else if (solver == "glpsol") {
            command << " --lp " << lpfile;
            if (timeout > 0) command << " --tmlim " << timeout;
            command << " -w " << solfile;
        } else if (solver == "gurobi_cl") { //Carmine 25.02.22 gurobi_cl included as MILP solver
            if (timeout > 0) command << " TimeLimit=" << timeout;
            command << " ResultFile=" << solfile;
            command << " " << lpfile;
        } else {
            assert(false);
        }
        return command.str();
    }

    static std::string createTempFilename(const std::string& prefix, const std::string& suffix = "") {
        char fname[1000];
        strcat(strcpy(fname, prefix.c_str()), ".XXXXXX");
        if (not suffix.empty()) strcat(fname, suffix.c_str());
        int file = mkstemps(fname, suffix.size());
        if (file == -1) return "";
        return fname;
    }

    static void deleteTempFilename(const std::string& filename) {
        unlink(filename.c_str());
    }
};
#endif // MILP_INTERFACE_H