// Copyright (C) 1996 DIMACS Center, Rutgers, The State University of New Jersey
// Author(s): John MacCuish (Los Alamos Nat. Lab.),

// This software is copyrighted by the DIMACS Center at Rutgers, The State
// University of New Jersey.  IT IS PROVIDED AS IS, AND THE AUTHORS, DIMACS, AND
// RUTGERS, THE STATE UNIVERSITY OF NEW JERSEY  DISCLAIM
// ALL LIABILITY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
// DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY
// DERIVATIVES THEREOF, EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

// THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,
// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE
// IS PROVIDED ON AN "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE
// NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
// MODIFICATIONS.

// The authors hereby grant permission to use, copy, modify, distribute,
// and license this software and its documentation for any purpose, provided
// that existing copyright notices are retained in all copies and that this
// notice is included verbatim in any distributions. No written agreement,
// license, or royalty fee is required for any of the authorized uses.
// Modifications to this software may be copyrighted by their authors
// and need not follow the licensing terms described here, provided that
// the new terms are clearly indicated on the first page of each file where
// they apply.

// Last File Update: 31-Jul-1996
// 


#include <math.h>
#include <stdio.h>
#include <LINK/graph/Graph.h>
#include <LINK/graph/Vertex.h>
#include <LINK/graph/Edge.h>
#include <LINK/graph/Attribute.h>
#include <LINK/basic/Queue.h>
#include <LINK/basic/Stack.h>
#include <LINK/basic/BinaryHeap.h>
#include <LINK/algorithm/Algorithms.h>

//
// get either edge flow or capacity assignments
//
void
edgeAssignments(Graph* graph, int** assignments, int assignmentflag)
{
    int i = 1;
    Edge *edge;
    const MSet<Edge*>& edges = graph->edges();
    Iterator<Edge*> get_edge(&edges);        
    while (get_edge(edge)) {
	Vertex* source = edge->sourceVertex();
	Vertex* sink = edge->sinkVertex();

	int j, k;
	getAttribute((GraphObject*) source, "index", j);
	getAttribute((GraphObject*) sink, "index", k);

	if (assignmentflag == 1)
	    getAttribute((GraphObject*) edge, "capacity", 
		assignments[j][k]);
	else {
	    if (assignments[j][k] < 0)
	        assignments[j][k] = 0;
	    cout << "Edge " << i <<
	        " has a flow value of " << assignments[j][k] << endl;
	    i++;
	}
    }
}


//
// initial exact labeling from sink
//
void
exactLabeling(Vertex** vertex_array, int sinkindex, int* height)
{
    // Exact labeling from sink
    Queue<Vertex*> BFSqueue;    
    Vertex *v;
    BFSqueue.enqueue(vertex_array[sinkindex]);

    Vertex* frontvertex;
    while (!BFSqueue.emptyQ()) {
	frontvertex = BFSqueue.front();
        Set<Vertex*> neighbors = frontvertex->inNeighbors();
        Iterator<Vertex*> get_vertices(&neighbors);
	Vertex* sinkvertex;	
	int frontindex;
	getAttribute((GraphObject*) frontvertex, "index", frontindex);
	
	while (get_vertices(sinkvertex)) {
	    int index;
	    getAttribute((GraphObject*) sinkvertex, "index", index);
	    if (height[index] == 0) {
		height[index] = height[frontindex] + 1;
		    BFSqueue.enqueue(sinkvertex);
	    }
	}
	v = BFSqueue.dequeue();
    }
}


//
// BFS from source or sink on residual network
//
void
sourceSinkBFS(Vertex** vertex_array, int sink_or_source, 
		int* height, int** residual_flow,int* heightcheck)
{
    Queue<Vertex*> BFSqueue;    
    Vertex* v;
    BFSqueue.enqueue(vertex_array[sink_or_source]);

    Vertex* frontvertex;
    while (!BFSqueue.emptyQ()) {
	frontvertex = BFSqueue.front();
        Set<Vertex*> neighbors = frontvertex->neighbors();
        Iterator<Vertex*> get_vertices(&neighbors);
	Vertex* sinkvertex;	
	int index_i;
	getAttribute((GraphObject*) frontvertex, "index", index_i);
	while (get_vertices(sinkvertex)) {
	    int index_j;
	    getAttribute((GraphObject*) sinkvertex, "index", index_j);
	    if (residual_flow[index_i][index_j] > 0 && heightcheck[index_i] == 0){
		height[index_i] = height[index_j] + 1;
		heightcheck[index_i] = 1;
		BFSqueue.enqueue(frontvertex);
	    }
	}
	v = BFSqueue.dequeue();
    }
}


void
globalRelabeling(int number_of_vertices, int sourceindex, 
		 int sinkindex, Vertex** vertex_array, int* height, 
		 int** residual_flow)
{
    int* heightcheck = new int[number_of_vertices];

    for (int i = 0; i < number_of_vertices; i++)
	heightcheck[i] = 0;
    heightcheck[sourceindex] = 1;
    heightcheck[sinkindex] = 1;

    sourceSinkBFS(vertex_array, sinkindex, height, residual_flow, heightcheck);
    sourceSinkBFS(vertex_array, sourceindex, height, residual_flow, heightcheck);
    delete [] heightcheck;
}
	

//
// push all flow from the source
//
void
pushSource(int use_key, Container<int>* container, int sinkindex,
	   int** flow, int** capacity, int* excess, Vertex* sourcevertex, 
	   int* datakey, int* height)
{
    Deque<int>* store;
    Dictionary<int,int>* dict;
    if (use_key)
	dict = (Dictionary<int,int>*) container;
    else
	store = (Deque<int>*) container;

    Set<Vertex*> neighbors = sourcevertex->outNeighbors();
    Iterator<Vertex*> get_vertices(&neighbors);
    Vertex* sinkvertex;
    while (get_vertices(sinkvertex)) {
	int j, k;
	getAttribute((GraphObject*) sourcevertex, "index", j);
	getAttribute((GraphObject*) sinkvertex, "index", k);
	flow[j][k] = excess[k] = capacity[j][k];
	flow[k][j] = -flow[j][k];
	if (k != sinkindex && excess[k] != 0)
	    if (use_key)
	        dict->insert(datakey[k], k);
	    else
	        store->insert(k);
	if (height[k] == 0)
	  height[k]++;
    }
}


void 
pushEdge(int* excess, int foundindex, int use_key, 
	      Container<int>* container,
	      Vertex** vertex_array, int** residual_flow, int* height,
	      int** flow,int sourceindex, int sinkindex, int &edgepushed,
	      int &maximum_distance, int* datakey, int &refpushcount)
{
    Deque<int>* store;
    Dictionary<int,int>* dict;
    if (use_key)
	dict = (Dictionary<int,int>*) container;
    else
	store = (Deque<int>*) container;

    Set<Vertex*> neighbors = vertex_array[foundindex]->neighbors();
    Iterator<Vertex*> get_vertex(&neighbors);
    Vertex* adjvertex;
    while (excess[foundindex] > 0 && (get_vertex(adjvertex))) {
	int i;
	getAttribute((GraphObject*) adjvertex, "index", i);
	if (residual_flow[foundindex][i] > 0) {
	    if (height[foundindex] == height[i] + 1) {
		int deltaflow;
		if (excess[foundindex] < residual_flow[foundindex][i])
		    deltaflow = excess[foundindex];
		else
		    deltaflow = residual_flow[foundindex][i];
		flow[foundindex][i] = flow[foundindex][i] + deltaflow;
		flow[i][foundindex] = -flow[foundindex][i];
		excess[foundindex] = excess[foundindex] - deltaflow;
		if (excess[foundindex] != 0)
		    if (use_key)
		        dict->insert(datakey[foundindex], foundindex);
		    else
		        store->insert(foundindex);
		excess[i] = excess[i] + deltaflow;
		if (i != sourceindex && i != sinkindex)
		    if (use_key)
		        dict->insert(datakey[i], i);
		    else
			store->insert(i);
		residual_flow[foundindex][i] =
		    residual_flow[foundindex][i] - deltaflow;
		residual_flow[i][foundindex] =
		    residual_flow[i][foundindex] + deltaflow;
		edgepushed++;
	       refpushcount++;
	    }
	    else if (edgepushed == 0 && maximum_distance > height[i])
		maximum_distance = height[i];
	}
    }
}


//
// push and relabel (lift)
//
void
pushRelabel(int use_key, Container<int>* container, Vertex** vertex_array, 
	    int sourceindex, int number_of_vertices, int sinkindex, 
	    int foundindex, int* height, int* excess, int** flow, 
	    int** residual_flow, int* datakey, int &refpushcount,int &refliftcount)
{
    while (excess[foundindex] > 0) {
	int edgepushed = 0;
	int maximum_distance = 2*number_of_vertices;
	pushEdge(excess, foundindex, use_key, container, vertex_array, 
		 residual_flow, height, flow, sourceindex, sinkindex,
		 edgepushed, maximum_distance, datakey, refpushcount);
	// lift (relabel) vertex
	if (edgepushed == 0){
	  refliftcount++;
	  height[foundindex] = maximum_distance + 1;
	}
    }
}


//
// push and relabel (lift) with gap relabeling
//
void
pushRelabelGap(int use_key, Container<int>* container, Vertex** vertex_array, 
	    int sourceindex, int number_of_vertices, int sinkindex, 
	    int foundindex, int* gapcount, int* height, int* excess, 
	    int** flow, int** residual_flow, 
	    int* heightcount, int* datakey, int &refpushcount,int &refliftcount)
{
    while (excess[foundindex] > 0) {
	int edgepushed = 0;
	int maximum_distance = 2*number_of_vertices;
	pushEdge(excess, foundindex, use_key, container, vertex_array, 
	 	 residual_flow, height, flow, sourceindex, sinkindex,
		 edgepushed, maximum_distance, datakey, refpushcount);
	// lift (relabel) vertex
	if (edgepushed == 0) {
	    // Check to see if a gap has been found
	    if (height[foundindex] < number_of_vertices &&
	       heightcount[height[foundindex]] == 1) {
		cout << "gap found" << endl;
		(*gapcount)++;
		int gap = height[foundindex];
		for (int i=0; i < number_of_vertices; i++) {
		    if (height[i] >= gap && height[i] < number_of_vertices) {
			heightcount[height[i]]--;
			height[i] = number_of_vertices;
			heightcount[height[i]]++;
		    }
		}
	    } else {
	        refliftcount++;
		heightcount[height[foundindex]]--;
		height[foundindex] = maximum_distance + 1;
		heightcount[height[foundindex]]++;
	    }
	}
    }
} 


void 
maxflowoutput(Graph* graph, int** flow, int* excess, int sinkindex,
		   int assignmentflag)
{
    cout << "The edge assignments are: " << endl;
    // set the assignment flag to print the edge flow assignments.
    assignmentflag = 0;
    edgeAssignments(graph, flow, assignmentflag);
    cout << "The maxflow is " << excess[sinkindex] << endl;
}

//***********************************************************************
//    vertexselection methods:
//    	LIFO (Stack)
//    	FIFO (Queue)
//    	MAXIMUM_EXCESS (Binary Heap)
//	HIGHEST_LABEL  (Binary Heap)
//
//    relabelchoice:
//	GLOBAL	: global relabeling
//	GAP	: gap relabeling
//	OTHER 	: neither
//
//    labelingchoice:
//	GLOBAL : global labeling
//	NAIVE  : naive
//***********************************************************************
    
int
GoldbergTarjan(Graph* graph,    int vertexselection=LIFO, 
				int relabelchoice=GLOBAL,
				int labelingchoice=GLOBAL)
{
    Deque<int>* store = 0;
    Dictionary<int,int>* dict = 0;

    int number_of_vertices = graph->order();
    int number_of_edges = graph->size();
    int minparameter = (int) sqrt(number_of_edges);
    int maxparameter = (number_of_edges/2) + 1;
    int discharge = 1;
    int globalparameter = 0;
    int use_key = 0;

    // initialize arrays and matrices
    int* height = new int[number_of_vertices];
    int* excess = new int[number_of_vertices];
    int** residual_flow = new int*[number_of_vertices];
    int** flow = new int*[number_of_vertices];
    int** capacity = new int*[number_of_vertices];
    int* datakey;

    switch (vertexselection) {
      case LIFO:
        store = new Stack<int>;
        datakey = excess;
        break;
      case FIFO:
        store = new Queue<int>;
        datakey = excess;
        break;
      case MAXIMUM_EXCESS:
        dict = new BinaryHeap<int, int>;
        datakey = excess;
	use_key = 1;
        break;
      default:
        dict = new BinaryHeap<int, int>;
        datakey = height;
	use_key = 1;
        break;
    }
    
    // Source vertex is the first vertex in vertex list and array
    const int sourceindex = 0;
    
    // Sink vertex is the last vertex in vertex list and array
    const int sinkindex = number_of_vertices-1; 

    // Create residual_flow, flow, and capacity adjacency matrices, initialize
    // height and excess arrays to 0
    int i;
    for (i=0; i < number_of_vertices; i++) {
        residual_flow[i] = new int[number_of_vertices];
	flow[i] = new int[number_of_vertices];
        capacity[i] = new int[number_of_vertices];
	height[i] = excess[i] = 0;
    }
    
    // if gap relabeling initialize the height count array
    int* heightcount = new int[number_of_vertices*2];
    if (relabelchoice == GAP) {
	for (i=0; i < number_of_vertices*2; i++)
	    heightcount[i] = 0;
    }

    // create vertex array
    Vertex** vertex_array; 

    // Initialize flow and capacity matrices to 0
    Attribute<int> index_attr(graph, "index", 0);
    for (i=0; i < number_of_vertices; i++) {
        vertex_array = graph->vertexStart();
	setAttribute((GraphObject*) vertex_array[i], "index", i);
	for (int j=0; j < number_of_vertices; j++)
	    flow[i][j] = capacity[i][j] = 0;
    }

    // Get edge capacities
    int assignmentflag = 1;
    edgeAssignments(graph, capacity, assignmentflag); 
    
    // Exact labeling from the sink to the source
    if (labelingchoice == GLOBAL)
        exactLabeling(vertex_array, sinkindex, height);

    // if gap relabeling update the height count array after the exactlabeling
    if (relabelchoice == GAP) {
	// height count array for gap relabeling
	for (i=0; i < number_of_vertices; i++)
	    heightcount[height[i]]++;
    }
    
    // raise the source to the number of vertices
    height[0] = number_of_vertices;

    // Push all flow from source along adjacent edges
    Vertex* sourcevertex = vertex_array[sourceindex];	
    if (use_key)
        pushSource(use_key, (Container<int>*) dict, 
	    sinkindex, flow, capacity, excess, sourcevertex, datakey, height);
    else
        pushSource(use_key, (Container<int>*) store, 
	    sinkindex, flow, capacity, excess, sourcevertex, datakey, height);

    // Initialize residual matrix after flow from source has been pushed
    for (i=0; i < number_of_vertices; i++)
	for (int j=0; j < number_of_vertices; j++)
	    residual_flow[i][j] = capacity[i][j] - flow[i][j];

    // Main loop of maxflow with three different choices of relabeling schemes
    int foundindex = 0;
    int liftcount = 0;
    int &refliftcount = liftcount;
    int pushcount = 0;
    int &refpushcount = pushcount;
    switch (relabelchoice) {
      // Global relabeling
      case GLOBAL:
	{
	cout << "Suggested range for global relabeling parameter: ";
	cout << minparameter << " to " << maxparameter << endl;
	cout << "Enter global parameter value:  ";
	cin >> globalparameter;
        int globalcount = 0;
	while (foundindex = store->get()) {
	    if (discharge%globalparameter == 0) {
		globalRelabeling(number_of_vertices, sourceindex, 
			sinkindex, vertex_array, height, residual_flow);
		cout << "global relabel" << endl;
		globalcount++;
	    }
	    if (use_key)
	        pushRelabel(use_key, (Container<int>*) dict, 
			vertex_array, sourceindex,
			number_of_vertices, sinkindex, foundindex,
			height, excess, flow, residual_flow,datakey, 
			refpushcount, refliftcount);
	    else
	        pushRelabel(use_key, (Container<int>*) store, 
			vertex_array, sourceindex,
			number_of_vertices, sinkindex, foundindex,
			height, excess, flow, residual_flow,datakey, 
			refpushcount, refliftcount);
	    discharge++;
	}
	maxflowoutput(graph, flow, excess, sinkindex, assignmentflag);
        cout << "Number of global relabels is:  " << globalcount << endl;
	cout << "Number of pushes is:  " << pushcount << endl;
        cout << "Number of lifts is:  " << liftcount << endl;
	break;
      }
      // Gap relabeling
      case GAP: {
	int gapcount = 0;
	while (foundindex = store->get()) {
	   if (use_key)
	       pushRelabelGap(use_key, (Container<int>*) dict, 
		   vertex_array, sourceindex,
		   number_of_vertices, sinkindex, foundindex, &gapcount,
		   height, excess, flow, residual_flow, heightcount, datakey,
		   refpushcount, refliftcount);
	    else
	        pushRelabelGap(use_key, (Container<int>*) store, 
		   vertex_array, sourceindex,
		   number_of_vertices, sinkindex, foundindex, &gapcount,
		   height, excess, flow, residual_flow, heightcount, datakey,
		   refpushcount, refliftcount);
	}
	maxflowoutput(graph, flow, excess, sinkindex, assignmentflag);
        cout << "Number of gap relabels is:  " << gapcount << endl;
	cout << "Number of pushes is:  " << pushcount << endl;
        cout << "Number of lifts is:  " << liftcount << endl;
	break;
      }
      // No relabeling scheme.
      default:
	while (foundindex = store->get()) {
	    if (use_key)
	        pushRelabel(use_key, (Container<int>*) dict, 
		    vertex_array, sourceindex, 
		    number_of_vertices, sinkindex, foundindex, height, excess, 
		    flow, residual_flow, datakey, refpushcount, refliftcount);
	    else
	        pushRelabel(use_key, (Container<int>*) store, 
		    vertex_array, sourceindex, 
		    number_of_vertices, sinkindex, foundindex, height, excess, 
		    flow, residual_flow, datakey, refpushcount, refliftcount);
	}
	maxflowoutput(graph, flow, excess, sinkindex, assignmentflag);
        cout << "Number of pushes is:  " << pushcount << endl;
        cout << "Number of lifts is:  " << liftcount << endl;
	break;
    }	

    Edge *edge;
    const MSet<Edge*>& edges = graph->edges();
    Iterator<Edge*> get_edge(&edges);

    while (get_edge(edge)) {
	Vertex* source = edge->sourceVertex();
	Vertex* sink = edge->sinkVertex();
	int j, k;
	getAttribute((GraphObject*) source, "index", j);
	getAttribute((GraphObject*) sink, "index", k);
	setAttribute((GraphObject*) edge,"flowvalue",flow[j][k]);
    }

    int max_flow = excess[sinkindex];

    delete [] height;
    delete [] excess;
    delete [] residual_flow;
    delete [] flow;
    delete [] capacity;
    delete [] heightcount;
    return max_flow;
}