#include "betweenness.h"

#include <boost/config.hpp>
#include <iostream>
#include <string>
#include <vector>
#include <stdio.h>
#include <boost/graph/edmunds_karp_max_flow.hpp> 
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_utility.hpp> 
#include <boost/graph/push_relabel_max_flow.hpp>
#include <boost/graph/read_dimacs.hpp>



int
Betweenness_Init (
    Tcl_Interp *interp) /* interpreter */
{


    Tcl_CreateCommand (interp, "bTwCompute", bTwComputeC, 
            (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);

    return TCL_OK;
}

/* Tcl command 
    bTwCompute adj_matrix */


int bTwComputeC (
    ClientData clientData,  /* NULL */
    Tcl_Interp *interp,     /* interpreter */
    int argc,               /* number of arguments */
    char *argv[])           /* argument list */
{ 
  
    int listArgc;
    char **listArgv;
    int numvert;
    long x;
    long maxFlow[256], vmaxFlow[256];

    using namespace boost;

    typedef adjacency_list_traits < vecS, vecS, directedS > Traits;
    typedef adjacency_list < vecS, vecS, directedS,
      property < vertex_name_t, std::string >,
      property < edge_capacity_t, long,
      property < edge_residual_capacity_t, long,
      property < edge_reverse_t, Traits::edge_descriptor > > > > Graph;


    if (argc != 2) {
	Tcl_AppendResult (interp, 
			  "wrong # args: should be \"",
			  argv[0], 
			  " adj_matrix\"",
			  (char *) NULL);
	return TCL_ERROR;
    }


    if (Tcl_SplitList (interp, argv[1], &listArgc, &listArgv) != TCL_OK) {
		Tcl_AppendResult (interp, 
			  "wrong first split # args: should be \"",
			  argv[0], 
			  " adj_matrix\"",
			  (char *) NULL);
	return TCL_ERROR;
    }

    Graph g;
    property_map < Graph, edge_capacity_t >::type
      capacity = get(edge_capacity, g);
    property_map < Graph, edge_reverse_t >::type 
	  rev = get(edge_reverse, g);         
	property_map < Graph, edge_residual_capacity_t >::type
      residual_capacity = get(edge_residual_capacity, g); 
	property_map<Graph, vertex_index_t>::type indexmap = get(vertex_index, g);

    typedef  graph_traits<Graph>::vertices_size_type vertices_size_type;
    typedef  graph_traits<Graph>::vertex_descriptor vertex_descriptor;
    typedef  graph_traits<Graph>::edge_descriptor edge_descriptor;

	edge_descriptor ed1, ed2;   
    std::vector<vertex_descriptor> verts;
 		  

    if (listArgc > 1) {
   
        bool inserted;
		

    	numvert = listArgc;
      
        for (long vi = 0; vi < numvert; vi++){
   	       verts.push_back(add_vertex(g));          
		}


    	for ( long i = 0; i < listArgc; i++) {
	        int list2Argc;
            char **list2Argv;

            if (Tcl_SplitList (interp, listArgv[i], &list2Argc, &list2Argv) 
		    != TCL_OK) {
                Tcl_Free ((char *) listArgv);
	        return TCL_ERROR;
            }
	        if (list2Argc != numvert) {
	            Tcl_AppendResult (interp, 
		     	  "\"", listArgv[i], "\"", " is not a symmetric matrix",
			      (char *) NULL);
                   Tcl_Free ((char *) list2Argv);
                   Tcl_Free ((char *) listArgv);
	               return TCL_ERROR;
			}

            long j;
            for (j = 0; j < list2Argc; j++) {
	           x = atoi (list2Argv[j]);
               if ((x > 0) && (i != j) ){

                   tie(ed1, inserted) = add_edge(verts[i], verts[j],g);
				   tie(ed2, inserted) = add_edge(verts[j], verts[i],g);

                   capacity[ed1] = x;
				   capacity[ed2] = 0;
                   rev[ed1] = ed2;
				   rev[ed2] = ed1;
				  
               }

            }
	   
            Tcl_Free ((char *) list2Argv);
        }
    }

    Tcl_Free ((char *) listArgv);

	 
    for(long m = 0; m < numvert; m++) {  
			   maxFlow[m] = 0;
               vmaxFlow[m] = 0;
	}

    for(long j = 0; j < numvert; j++) {

        for (long k = 0; k < numvert; k++) {

			if((j != k)) { 
				              
/*			
			   std::vector<default_color_type> color(num_vertices(g));
               std::vector<Traits::edge_descriptor> pred(num_vertices(g));
               long flow = edmunds_karp_max_flow
               (g, j, k, capacity, residual_capacity, rev, &color[0], &pred[0]);
*/

               property_map<Graph, vertex_index_t>::type indexmap = get(vertex_index, g);
               long flow = push_relabel_max_flow(g, verts[j], verts[k],capacity, residual_capacity, rev,
			   	                                indexmap);

    			  
               for(int i = 0; i < numvert; i++) {
                 if(!(( k ==i) || (j == i))) {

                    maxFlow[i] = maxFlow[i] + flow;
    
                    graph_traits < Graph >::out_edge_iterator ei, e_end;                   
                    for (tie(ei, e_end) = out_edges(verts[i], g); ei != e_end; ++ei)      
				       if (capacity[*ei] > 0) {
                          long btw;                     
				   	      btw = (capacity[*ei] - residual_capacity[*ei]);
                          if (btw > 0) vmaxFlow[i] = vmaxFlow[i] +  btw;
					   }

					}
			   }  /* i */ 
			} /*j != k */

        } /* k */

    } /* j */

    for(long i = 0; i < numvert; i++) {  
      char buf[256];
      sprintf (buf, "{%d %d %f} ", 
			      maxFlow[i], vmaxFlow[i], (float)vmaxFlow[i]/maxFlow[i]);
      Tcl_AppendResult (interp, buf, (char *) NULL);
	}

    return TCL_OK;

}