00001 #include<scil/scil.h>
00002 #include<scil/constraints/SteinerArborescence.h>
00003 #include <boost/graph/adjacency_list.hpp>
00004 #include <boost/graph/iteration_macros.hpp>
00005
00006 using namespace std;
00007 using namespace SCIL;
00008 using namespace boost;
00009
00010
00011 typedef boost::adjacency_list<vecS,vecS,bidirectionalS> Graph;
00012 typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
00013 typedef boost::graph_traits<Graph>::edge_descriptor edge_descriptor;
00014
00015
00016 void read_network(char* file_name, Graph& G, map<edge_descriptor, double>& edge_weights, map<vertex_descriptor, pair<int, int> >& vertexmap)
00017 {
00018 ifstream file;
00019 file.open(file_name);
00020
00021 int number_of_vertices;
00022
00023 file >> number_of_vertices;
00024 double x,y;
00025 vertex_descriptor v,w;
00026
00027 for(int i = 0; i<number_of_vertices; i++){
00028 file >> x >> y;
00029 vertexmap[add_vertex(G)] = make_pair(x,y);
00030 }
00031
00032 file.close();
00033
00034 BGL_FORALL_VERTICES(v, G, Graph){
00035 BGL_FORALL_VERTICES(w, G, Graph){
00036 if(v != w){
00037 edge_descriptor e = add_edge(v,w,G).first;
00038 double dist = sqrt(((vertexmap[w].second - vertexmap[v].second) * (vertexmap[w].second - vertexmap[v].second) ) + ((vertexmap[w].first - vertexmap[v].first) * (vertexmap[w].first - vertexmap[v].first)));
00039 edge_weights[e] = dist;
00040 }
00041 }
00042 }
00043 }
00044
00045
00046 int main(int argc, char** argv){
00047
00048 Graph G;
00049 map<edge_descriptor, double> edge_weights;
00050 map<vertex_descriptor, pair<int, int> > vertexmap;
00051 var_map<edge_descriptor> VM;
00052
00053 read_network(argv[1], G, edge_weights, vertexmap);
00054
00055 vertex_descriptor root = *(vertices(G).first);
00056 list<vertex_descriptor> terminals;
00057 int num_terminals = atoi(argv[2]);
00058
00059 int i = 0;
00060 BGL_FORALL_VERTICES(v, G, Graph){
00061 if(v != root && (i < num_terminals || num_terminals == 0)){
00062 terminals.push_back(v);
00063 i++;
00064 }
00065 }
00066
00067 ILP_Problem IP(Optsense_Min);
00068
00069 BGL_FORALL_EDGES(e, G, Graph)
00070 VM[e] = IP.add_binary_variable(edge_weights[e]);
00071
00072 IP.add_sym_constraint(new SteinerArborescence<Graph>(G, root, terminals, VM));
00073
00074 IP.optimize();
00075
00076 cout<<"Optimum: "<<IP.get_optimum()<<endl;
00077 BGL_FORALL_EDGES(e, G, Graph){
00078 if(IP.get_solution(VM[e]) == 1)
00079 cout<<e<<endl;
00080 }
00081 }
00082
