miniksp.cc

   1
   2 ///////////////////////////////////////////////////////////////////////////
   3 // START_IP_HEADER                                                       //
   4 //                                                                       //
   5 // This program is free software: you can redistribute it and/or modify  //
   6 // it under the terms of the version 3 of the GNU General Public License //
   7 // as published by the Free Software Foundation.                         //
   8 //                                                                       //
   9 // This program is distributed in the hope that it will be useful, but   //
  10 // WITHOUT ANY WARRANTY; without even the implied warranty of            //
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU      //
  12 // General Public License for more details.                              //
  13 //                                                                       //
  14 // You should have received a copy of the GNU General Public License     //
  15 // along with this program. If not, see <http://www.gnu.org/licenses/>.  //
  16 //                                                                       //
  17 // Written by and Copyright (C) Francois Fleuret                         //
  18 // Contact <francois.fleuret@idiap.ch> for comments & bug reports        //
  19 //                                                                       //
  20 // END_IP_HEADER                                                         //
  21 ///////////////////////////////////////////////////////////////////////////
  22
  23 #define VERBOSE
  24
  25 #include <iostream>
  26 #include <fstream>
  27 #include <cmath>
  28 #include <stdio.h>
  29 #include <stdlib.h>
  30 #include <float.h>
  31
  32 using namespace std;
  33
  34 typedef float scalar_t;
  35
  36 #ifdef DEBUG
  37 #define ASSERT(x) if(!(x)) { \
  38   std::cerr << "ASSERT FAILED IN " << __FILE__ << ":" << __LINE__ << endl; \
  39   abort(); \
  40 }
  41 #else
  42 #define ASSERT(x)
  43 #endif
  44
  45 void raise_es(int nb_edges, scalar_t *el) {
  46   scalar_t min_es = el[0];
  47   for(int e = 1; e < nb_edges; e++) {
  48     min_es = min(min_es, el[e]);
  49   }
  50   for(int e = 0; e < nb_edges; e++) {
  51     el[e] -= min_es;
  52   }
  53 }
  54
  55 void add_dpsi_es(int nb_edges, scalar_t *el, int *ea, int *eb, scalar_t *psi) {
  56   for(int e = 0; e < nb_edges; e++) {
  57     el[e] += psi[ea[e]] - psi[eb[e]];
  58   }
  59 }
  60
  61 void find_shortest(int nb_vertices,
  62                    int nb_edges, scalar_t *el, int *ea, int *eb,
  63                    int source, int sink,
  64                    int *edge_back, scalar_t *dist) {
  65   for(int v = 0; v < nb_vertices; v++) {
  66     // dist[v] = FLT_MAX;
  67     dist[v] = 999;
  68   }
  69
  70   dist[source] = 0;
  71
  72   for(int v = 0; v < nb_vertices; v++) {
  73     edge_back[v] = -1;
  74   }
  75
  76   for(int e = 0; e < nb_edges; e++) {
  77     ASSERT(el[e] >= 0);
  78   }
  79
  80   int nb_changes;
  81   scalar_t d;
  82   do {
  83     nb_changes = 0;
  84     for(int e = 0; e < nb_edges; e++) {
  85       d = dist[ea[e]] + el[e];
  86       if(d < dist[eb[e]]) {
  87         nb_changes++;
  88         dist[eb[e]] = d;
  89         edge_back[eb[e]] = e;
  90       }
  91     }
  92   } while(nb_changes > 0);
  93 }
  94
  95 void find_best_paths(int nb_vertices,
  96                      int nb_edges, scalar_t *el, int *ea, int *eb,
  97                      int source, int sink,
  98                      int *edge_occupation) {
  99   scalar_t *dist = new scalar_t[nb_vertices];
 100   int *edge_back = new int[nb_vertices];
 101   scalar_t *tmp_el = new scalar_t[nb_edges];
 102
 103   for(int e = 0; e < nb_edges; e++) {
 104     tmp_el[e] = el[e];
 105     edge_occupation[e] = 0;
 106   }
 107
 108   raise_es(nb_edges, tmp_el);
 109
 110   scalar_t s;
 111   do {
 112     find_shortest(nb_vertices, nb_edges, tmp_el, ea, eb, source, sink, edge_back, dist);
 113     add_dpsi_es(nb_edges, tmp_el, ea, eb, dist);
 114     s = 0.0;
 115     if(edge_back[sink] >= 0) {
 116       int v = sink;
 117       while(v != source) {
 118         int e = edge_back[v];
 119         ASSERT(eb[e] == v);
 120         v = ea[e];
 121 #ifdef VERBOSE
 122         cout << "INVERTING " << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
 123 #endif
 124         s += el[e];
 125         int k = eb[e];
 126         eb[e] = ea[e];
 127         ea[e] = k;
 128         tmp_el[e] = - tmp_el[e];
 129         edge_occupation[e] = 1 - edge_occupation[e];
 130       }
 131 #ifdef VERBOSE
 132       cout << "FOUND A PATH OF LENGTH " << s << endl;
 133 #endif
 134     }
 135   } while(s < 0);
 136
 137   delete[] tmp_el;
 138   delete[] dist;
 139   delete[] edge_back;
 140 }
 141
 142 int main(int argc, char **argv) {
 143
 144   if(argc < 2) {
 145     cerr << argv[0] << " <graph file>" << endl;
 146     exit(EXIT_FAILURE);
 147   }
 148
 149   ifstream file(argv[1]);
 150
 151   int nb_edges, nb_vertices;
 152   int source, sink;
 153
 154   if(file.good()) {
 155     file >> nb_vertices >> nb_edges;
 156     file >> source >> sink;
 157
 158     cout << "nb_edges = " << nb_edges << endl;
 159     cout << "nb_vertices = " << nb_vertices << endl;
 160     cout << "source = " << source << endl;
 161     cout << "sink = " << sink << endl;
 162
 163     scalar_t *el = new scalar_t[nb_edges];
 164     int *ea = new int[nb_edges];
 165     int *eb = new int[nb_edges];
 166     int *edge_occupation = new int[nb_edges];
 167
 168     for(int e = 0; e < nb_edges; e++) {
 169       file >> ea[e] >> eb[e] >> el[e];
 170 #ifdef VERBOSE
 171       cout << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
 172 #endif
 173     }
 174
 175     find_best_paths(nb_vertices, nb_edges, el, ea, eb, source, sink,
 176                     edge_occupation);
 177
 178     for(int e = 0; e < nb_edges; e++) {
 179       if(edge_occupation[e]) {
 180         cout << ea[e] << " " << eb[e] << endl;
 181       }
 182     }
 183
 184     delete[] edge_occupation;
 185     delete[] el;
 186     delete[] ea;
 187     delete[] eb;
 188
 189   } else {
 190
 191     cerr << "Can not open " << argv[1] << endl;
 192     exit(EXIT_FAILURE);
 193
 194   }
 195
 196   exit(EXIT_SUCCESS);
 197 }