working on compact MIP

examples/Optimization/MINLP/LABS/LABS.cpp

@@ -16,13 +16,27 @@ using namespace gravity;
 using namespace std;
 
 int main(int argc, char * argv[]){
+
+    double total_time_start = get_wall_time();
+
     /* binary string dimension*/
     int n = 3;
+    bool new_algorithm = false;
+    bool nlp = false;
+    string alg = "mip";
     if(argc>=2){
         n=stoi(argv[1]);
     }
+    if(argc>=3){
+        alg = argv[2];
+    }
+
+    if(alg=="new")
+        new_algorithm =true;
+    if(alg=="nlp")
+        nlp =true;
     DebugOn("Optimizing for N = " << n << endl);
-    bool new_algorithm = false;
+
     if(new_algorithm){
         vector<Model<>> models;
         for (int i = 0; i<n-1; i++) {
@@ -44,33 +58,84 @@ int main(int argc, char * argv[]){
     var<> c("c");
     indices s_ids = range(0,n-1);
     indices c_ids = range(1,n-1);
-    bool unconstrained = true;
-    if(unconstrained){
+    int opt_obj = 0;
+    if(nlp){
         s.exclude_zero();
         M_obj.add(s.in(s_ids));
         M_obj.add(c.in(c_ids));
+        indices pairs("pairs"), quad_terms("quad_terms"), multi_terms("multi_terms"), multi_quad_terms("multi_quad_terms");
         func<int> obj;
         for (int k = 1; k<=n-1; k++) {
             func<int> cterm;
             for (int i = 0; i<n-k; i++) {
                 cterm += (s(to_string(i)))*(s(to_string(i+k)));
             }
-            cterm.eval_all();
-            Constraint<> C_sq_k("C_sq_"+to_string(k));
-            C_sq_k = c(k) - pow(cterm,2);
-            M_obj.add(C_sq_k>=0);
-//            obj += pow(cterm,2);
+            obj += pow(cterm,2);
         }
 //        obj.print();
-
-        M_obj.min(sum(c));
+        int nb_quad = obj._qterms->size();
+        DebugOn("Number of quadratic terms = " << nb_quad << endl);
+        int nb_mult = obj._pterms->size();
+        DebugOn("Number of multilinear terms = " << nb_mult << endl);
+
+        for(auto p: *obj._qterms)
+        {
+            auto idx1 = p.second._p->first->_indices->_ids->front().at(0);
+            auto idx2 = p.second._p->second->_indices->_ids->front().at(0);
+            pairs.add(to_string(idx1)+","+to_string(idx2));
+        }
+        quad_terms.is_subset(pairs);
+        multi_quad_terms.is_subset(pairs);
+        string pair_idx;
+        int nb_row = 0;
+        for(auto p: *obj._pterms)
+        {
+            set<int> tabu;
+            for(auto it = p.second._l->begin(); it != p.second._l->end(); it++) {
+                bool existing_pair = false;
+                auto idx1 = it->first->_indices->_ids->front().at(0);
+                if(tabu.count(idx1)>0)
+                    continue;
+                for(auto it2 = next(it); it2 != p.second._l->end(); it2++) {
+                    auto idx2 = it2->first->_indices->_ids->front().at(0);
+                    if(tabu.count(idx2)>0)
+                        continue;
+                    pair_idx = to_string(idx1)+","+to_string(idx2);
+                    if(pairs.has_key(pair_idx)){
+                        existing_pair = true;
+                        multi_quad_terms.add_in_row(nb_row, pair_idx);
+                        tabu.insert(idx1);
+                        tabu.insert(idx2);
+                        break;
+                    }
+                }
+                if(existing_pair){
+                    continue;
+                }
+                else{
+                    pairs.add(pair_idx);
+                }
+            }
+            nb_row++;
+        }
+        DebugOn("Number of quadratic terms after decomposing quartic expressions = " << pairs.size() << endl);
+        pairs.print();
+
+        M_obj.min(obj);
         M_obj.print();
+//        auto g = M_obj.get_interaction_graph();
+//        g.print();
+//        g.get_tree_decomp_bags();
+//        auto ConvM = M_obj.relax();
+//        ConvM->print();
         s.initialize_binary();
+//        s.initialize_all(1);
         solver<> g_sol(M_obj,ipopt);
         g_sol.run();
         M_obj.print_solution();
         M_obj.round_solution();
         M_obj.print_solution();
+        opt_obj = M_obj.get_obj_val();
     }
     else{
         indices z_ids("z_ids");
@@ -122,6 +187,11 @@ int main(int argc, char * argv[]){
 //        M.print_solution();
 //        M.round_solution();
         M.print_solution();
+        opt_obj = M.get_obj_val();
     }
+    double total_time_end = get_wall_time();
+    auto total_time = total_time_end - total_time_start;
 
+    DebugOn("Optimal objective for n = " << n << " is " << opt_obj << endl);
+    DebugOn("Wall clock time = " << total_time << endl);
 }

src/GurobiProgram.cpp

@@ -503,7 +503,7 @@ bool GurobiProgram::solve(bool relax, double mipgap, double time_limit){
     //    grb_mod->getEnv().set(GRB_IntParam_NodeMethod, 1);
     grb_mod->getEnv().set(GRB_IntParam_LazyConstraints, 1);
 //    grb_mod->set(GRB_IntParam_Threads, 8);
-//    grb_mod->set(GRB_IntParam_NonConvex,2);
+    grb_mod->set(GRB_IntParam_NonConvex,2);
 //    grb_mod->set(GRB_DoubleParam_IntFeasTol, 1e-9);
     //       grb_mod->set(GRB_IntParam_NumericFocus,3);
     //     grb_mod->set(GRB_IntParam_PreCrush,0);

src/Net.cpp

@@ -545,16 +545,16 @@ void Net::get_tree_decomp_bags() {
         Debug(graph_clone->nodes.size() << endl);
         pair<string,vector<Node*>> bag_copy;
         pair<string,vector<Node*>> bag;
-        DebugOff("new bag = { ");
+        DebugOn("new bag = { ");
         for (auto nn: n->get_neighbours()) {
             if(!nn.second->_active) continue;
             bag_copy.second.push_back(nn.second);
             bag_copy.first += "," +nn.first;
             bag.second.push_back(get_node(nn.second->_name)); // Note it takes original node.
             bag.first += "," +nn.first;
-            DebugOff(nn.second->_name << ", ");
+            DebugOn(nn.second->_name << ", ");
         }
-        DebugOff(n->_name << "}\n");
+        DebugOn(n->_name << "}\n");
         graph_clone->remove_end_node();
         bag_copy.second.push_back(n);
         bag_copy.first += "," +n->_name;
@@ -610,9 +610,9 @@ void Net::get_tree_decomp_bags() {
     }
 
     delete graph_clone;
-    if(!sdp_3d_cuts){
+//    if(!sdp_3d_cuts){
         remove_redundant_bags();
-    }
+//    }
 
 }