Quadratic_Minimum_Spanning_Tree

This example shows how to compute a minimum spanning tree in a given graph G using the symbolic constraint SCIL::SpanTree.

In this case we also add quadratic terms to the objective function.

We first create a new instance IP of an SCIL::ILP_Problem as a minimization problem.

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

We then create an empty graph and initialize it using the function read_inst. This function also adds all required quadratic terms using IP.add_polynomial . The var_map VM is used to store the variables that are added to the model. These variables are binary and correspond to the edges of the graph G.

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

The edges corresponding to variables with value 1 in an optimal solution should form a valid spanning tree in G. To achieve this we add the symbolic constraint SCIL::SpanTree to the model.

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

In this example we also use a primal_heuristic. A primal heuristic tries to generate a primal solution from the fractional lp solution. This is usefull to speed up the branch&cut process. We add the previously implemented primal heuristc qmst_heuristic using IP.set_primal_heuristic.

   IP.set_primal_heuristic(new qmst_heuristic(G, VM));

Our model is now complete and we can compute an optimum solution by calling optimize.

   IP.optimize();

The optimal solution value can be accessed directedly using SCIL::ILP_Problem::get_optimum().

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

Finally we output the edges, which form a minimum spanning tree, using IP.get_solution().

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

The primal heuristic qmst_heuristic

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;
      }
      
};