qmst.cc

#include<iostream>
#include<scil/scil.h>
#include<scil/constraints/spantree.h>
#include<scil/core/monomial.h>
#include <boost/foreach.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/graph/kruskal_min_spanning_tree.hpp>

using namespace SCIL;
using namespace std;
using namespace boost;

typedef adjacency_list<vecS, vecS, undirectedS> Graph;
typedef boost::graph_traits<Graph> GraphTraits;
typedef GraphTraits::vertex_descriptor vertex_descriptor;
typedef GraphTraits::edge_descriptor edge_descriptor;

class qmst_heuristic : public primal_heuristic {
   private:
      Graph& G;
      var_map<edge_descriptor>& VM;

   public:
      qmst_heuristic(Graph& G_, var_map<edge_descriptor>& VM_)
         : G(G_), VM(VM_)
      { }

      int heuristic(subproblem& S, solution& newSol) {
         edge_descriptor e;
         map<edge_descriptor, int> tc;
         BGL_FORALL_EDGES(e, G, Graph) {
            tc[e] = int(10. * (1. - (S.value(VM[e]))));
         }

         list<edge_descriptor> tree;

         typedef map<edge_descriptor, int> weight_map_t;
         typedef map<vertex_descriptor, size_t> vertex_index_t;   

         vertex_index_t vertex_index_map;
         int i = 0;
         BGL_FORALL_VERTICES(v, G, Graph)
            vertex_index_map[v] = i++;
         
         associative_property_map<weight_map_t> pm_weight(tc);
         associative_property_map<vertex_index_t> pm_vertex_index(vertex_index_map);

         kruskal_minimum_spanning_tree(G, back_inserter(tree), boost::weight_map(pm_weight).vertex_index_map(pm_vertex_index));

         BGL_FORALL_EDGES(e, G, Graph) {
            newSol.set_value(VM[e], 0.);
         }
         BOOST_FOREACH(e, tree) {
            newSol.set_value(VM[e], 1.);
         }
         S.GM.extendSolution(newSol);

         return 1;
      }
      
};

void read_inst(const char* file_name, Graph& G, var_map<edge_descriptor>& VM, ILP_Problem& IP) {

   ifstream file;
   file.open(file_name);

   int m, n, h1, h2, t1, t2, w;
   file >> n >> m;
   vector<vertex_descriptor> nm(n+1);
   edge_descriptor* em = new edge_descriptor[n*n];

   //initialize node array
   for( int i = 1; i <= n; i++ )
      nm[i] = add_vertex(G);

   //read edges from file and add a corresponding variable to the IP
   for( int i = 0; i < m; i++ ) {
      file >> t1 >> h1 >> w;
      em[(t1-1)*n+(h1-1)] = add_edge(nm[t1], nm[h1], G).first;
      VM[em[(t1-1)*n+(h1-1)]] = IP.add_binary_variable(w);
   }

   //read values of the nonlinear part of objective function and add corresponding monomials to the IP
   for( int i = 0; i < m*(m-1); i++ ) {
      file >> t1 >> h1 >> t2 >> h2 >> w;
      IP.add_polynomial(VM[em[(t1-1)*n+(h1-1)]] *  VM[em[(t2-1)*n+(h2-1)]] *  w);
   }

   file.close();
   return;
}

int main(int argc, char** argv) 
{
   ILP_Problem IP(Optsense_Min);

   Graph G;
   var_map<edge_descriptor> VM;
   read_inst(argv[1], G, VM, IP);

   
   IP.add_sym_constraint(new SpanTree<Graph>(G, VM));
   
   IP.set_primal_heuristic(new qmst_heuristic(G, VM));

   IP.optimize();

   cout<<"optimal solution: "<<IP.get_optimum()<<endl;

   cout<<"spanning tree: "<<endl;
   BGL_FORALL_EDGES(e, G, Graph)
      if(IP.get_solution(VM[e]))
         cout<<e<<endl;
   
   return 0;
}

---