#ifndef __VBD__
#define __VBD__

#include <iostream>
#include <string>
#include <vector>
#include <map>
#include <set>
#include <cstdio>
#include <ctime>
#include <list>
#include <algorithm>    // std::random_shuffle
#include <queue>
#include <sstream>
#include <cctype>
#include "helper.h"
#include "Graph.h"
#include "makeGrid.h"
#include "VolcanoBlackholeDetection.h"
using namespace std;

inline void appendToFile(const char* filename, vector<int> &blackhole, Graph &g, int k, int i, int j, int tempBlackholesFlow, vector<int> &inFlow, vector<int> &outFlow)
{
	FILE* file = fopen(filename, "a");
	fprintf(file, "command=%d day=%d time frame=%d flow=%d size=%d\n", k, i, j, tempBlackholesFlow, blackhole.size());
	for(int i = 0; i < blackhole.size(); ++ i) {
		int nodef = g.eg[blackhole[i]].nodef;
		int nodet = g.eg[blackhole[i]].nodet;
		double nodef_lat = g.sg[nodef].la;
		double nodef_long = g.sg[nodef].lo;
		double nodet_lat = g.sg[nodet].la;
		double nodet_long = g.sg[nodet].lo;
		fprintf(file, "%d %d %d %d %lf %lf %d %lf %lf\n", i, blackhole[i], inFlow[blackhole[i]] - outFlow[blackhole[i]], nodef, nodef_lat, nodef_long, nodet, nodet_lat, nodet_long);
	}

	fprintf(file,"\n");
	fclose(file);
}

//actually append
inline void printVolcanoToFile(vector<int> &volcano, Graph &g, int k, int i, int j, int tempVolcanosFlow, vector<int> &inFlow, vector<int> &outFlow, const char* filename)
{
	//in/out is reversed for volcano use
	appendToFile(filename, volcano, g, k, i, j, tempVolcanosFlow, outFlow, inFlow);
}

//actually append
void printBlackholeToFile(vector<int> &blackhole, Graph &g, int k, int i, int j, int tempBlackholesFlow, vector<int> &inFlow, vector<int> &outFlow, const char* filename)
{
	appendToFile(filename, blackhole, g, k, i, j, tempBlackholesFlow, inFlow, outFlow);
}


struct data
{
	int edgeID;
	int absFlow;
};

struct LengthFlowRatio
{
	double length;
	int absflow;
	double ratio;
};

// random generator function:
inline int myrandom (int n)
{ 
	return myrand()%n;
}

bool checkAndUpdateArea(double &minLat, double &minLong, double &maxLat, double &maxLong, double fLat, double fLong, double tLat, double tLong, double areaThreshold)
{
	double oldminLat = minLat, oldminLong = minLong, oldmaxLat = maxLat, oldmaxLong = maxLong;

	updateInterval(minLat, maxLat, fLat);
	updateInterval(minLong, maxLong, fLong);
	updateInterval(minLat, maxLat, tLat);
	updateInterval(minLong, maxLong, tLong);

	double row = Geodist(minLat, minLong, minLat, maxLong) / 1000.0;
	double col = Geodist(maxLat, minLong, maxLat, maxLong) / 1000.0;
	double diagnose = Geodist(maxLat, maxLong, minLat, minLong) / 1000.0;
	double area = row * col;
	double diagnoseThreshold = sqrt(2.0 * areaThreshold);
	//printf("area=%lf areaThreshold=%lf diagnose=%lf diagnoseThreshold=%lf\n",area,areaThreshold,diagnose, diagnoseThreshold);

	if(area > areaThreshold || diagnose > diagnoseThreshold) {
		//restore
		minLat = oldminLat;
		minLong = oldminLong;
		maxLat = oldmaxLat;
		maxLong = oldmaxLong;
		return false;
	} else {
		return true;
	}
}

inline bool compareByDescAbsFlow(const data &a, const data &b)
{
	return a.absFlow > b.absFlow;
}

void UpdateBlackholeEdge(vector<int> &blackhole, vector<int> &blackholeEdge)
{
	for(int i = 0; i < blackhole.size(); ++ i) {
		blackholeEdge[blackhole[i]] = 1;
	}
}

void RemoveBlackholeEdgefromEdgeChoice2(vector<int> &blackhole, vector<int> &edgeChoice)
{
	for(int i = 0; i < blackhole.size(); ++ i) {
		edgeChoice[blackhole[i]]=1;
	}
}

// make sure (i, j) is not on the boundary!!
int getBlackholeUpperBound(vector<vector<gridCell>> &Grid, int i, int j)
{
	const int dx[8] = {-1, -1, 0, 1, 1, 1, 0, -1};
	const int dy[8] = {0, 1, 1, 1, 0, -1, -1, -1};
	int maximum = 0;
	for (int startDirection = 0; startDirection < 8; startDirection += 2) {
		int sum = 0;
		for (int iter = 0; iter < 3; ++ iter) {
			sum += Grid[i + dx[startDirection + iter & 7]][j + dy[startDirection + iter & 7]].BlackholeresidualFlow;
		}
		maximum = max(maximum, sum);
	}
	maximum += Grid[i][j].BlackholeresidualFlow;
	return maximum;
}

// make sure (i, j) is not on the boundary!!
int getVolcanoUpperBound(vector<vector<gridCell>> &Grid, int i, int j)
{
	const int dx[8] = {-1, -1, 0, 1, 1, 1, 0, -1};
	const int dy[8] = {0, 1, 1, 1, 0, -1, -1, -1};
	int maximum = 0;
	for (int startDirection = 0; startDirection < 8; startDirection += 2) {
		int sum = 0;
		for (int iter = 0; iter < 3; ++ iter) {
			sum += Grid[i + dx[startDirection + iter & 7]][j + dy[startDirection + iter & 7]].VolcanoresidualFlow;
		}
		maximum = max(maximum, sum);
	}
	maximum += Grid[i][j].VolcanoresidualFlow;
	return maximum;
}

//this one not use map (should be faster)
void blackholedynamicallyComputeUpperBound(Graph &g, vector<int> &inFlow, vector<int> &outFlow, const int flowThreshold, const double areaThreshold, const int numEdge ,
											vector<vector<gridCell>> &Grid, vector<int> &pickedEdge, vector<int> &edgeChoice, vector<int> &Blackhole,
											double maxLat, double maxLong, double minLat, double minLong)
{
	if (Grid.size() == 0) {
		debug("Grid is empty!!!");
		return;
	}
	int NumGridCellcol = Grid[0].size(), NumGridCellrow = Grid.size();
#define GRID_ID(i,j) ((i) * NumGridCellcol + (j))
    GridCellDistributor gridCellDistributor(minLat, maxLat, minLong, maxLong, NumGridCellcol, NumGridCellrow);
	int prunedCount=0;
	vector<int> affectedGridCell;
	for(int i = 0; i < Blackhole.size(); ++ i) {
		int row1, col1, row2, col2;
		int Nodef = g.eg[Blackhole[i]].nodef, Nodet = g.eg[Blackhole[i]].nodet;

		gridCellDistributor.findGridCell(row1, col1, g.sg[Nodef].la, g.sg[Nodef].lo);
		gridCellDistributor.findGridCell(row2, col2, g.sg[Nodet].la, g.sg[Nodet].lo);

		int delta = max(inFlow[Blackhole[i]] - outFlow[Blackhole[i]], 0);
		Grid[row1][col1].BlackholeresidualFlow -= delta;
		affectedGridCell.push_back(GRID_ID(row1, col1));
		if (row2 != row1 || col2 != col1) {
			Grid[row2][col2].BlackholeresidualFlow -= delta;
			affectedGridCell.push_back(GRID_ID(row2, col2));
		}
	}
#undef GRID_ID(i,j)

	sort(affectedGridCell.begin(), affectedGridCell.end());
	affectedGridCell.erase(unique(affectedGridCell.begin(), affectedGridCell.end()), affectedGridCell.end());

	//if max of sum of residual flow of 4 adj. grid cell < flow threshold, remove all edges of thegrid cell in edge choice
	for(int k = 0; k < affectedGridCell.size(); ++ k) {
		int cellId = affectedGridCell[k];
		int i = cellId / NumGridCellcol;
		int j = cellId % NumGridCellcol;

		if (i == 0 || j == 0 || i + 1 == NumGridCellrow || j + 1 == NumGridCellcol) {
			//it is on the boundary
			continue;
		}
		
		if (getBlackholeUpperBound(Grid, i, j) < flowThreshold) {
			Grid[i][j].noBlackhole = 1;
			if(Grid[i][j].noBlackhole == 1 && Grid[i][j].noVolcano == 1) {
				for(map<int, int>::iterator it = Grid[i][j]._Edge.begin(); it != Grid[i][j]._Edge.end(); ++ it) {
					prunedCount++;
					edgeChoice[(*it).first]=1;
				}
			}
		}
	}
	fprintf(stderr, "Blackhole: Dynamic Upper Bound:: pruned=%d\n",prunedCount);
}

//this one not use map (should be faster)
void volcanodynamicallyComputeUpperBound(Graph &g, vector<int> &inFlow, vector<int> &outFlow, const int flowThreshold, const double areaThreshold, const int numEdge,
											vector<vector<gridCell>> &Grid, vector<int> &pickedEdge, vector<int> &edgeChoice, vector<int> &Blackhole,
											double maxLat, double maxLong, double minLat, double minLong)
{
	if (Grid.size() == 0) {
		debug("Grid is empty!!!");
		return;
	}
	int NumGridCellcol = Grid[0].size(), NumGridCellrow = Grid.size();
#define GRID_ID(i,j) ((i) * NumGridCellcol + (j))
    GridCellDistributor gridCellDistributor(minLat, maxLat, minLong, maxLong, NumGridCellcol, NumGridCellrow);
	int prunedCount=0;
	vector<int> affectedGridCell;
	for(int i = 0; i < Blackhole.size(); ++ i) {
		int row1, col1, row2, col2;
		int Nodef = g.eg[Blackhole[i]].nodef, Nodet = g.eg[Blackhole[i]].nodet;

		gridCellDistributor.findGridCell(row1, col1, g.sg[Nodef].la, g.sg[Nodef].lo);
		gridCellDistributor.findGridCell(row2, col2, g.sg[Nodet].la, g.sg[Nodet].lo);

		int delta = max(outFlow[Blackhole[i]] - inFlow[Blackhole[i]], 0);
		Grid[row1][col1].VolcanoresidualFlow -= delta;
		affectedGridCell.push_back(GRID_ID(row1, col1));
		if (row2 != row1 || col2 != col1) {
			Grid[row2][col2].VolcanoresidualFlow -= delta;
			affectedGridCell.push_back(GRID_ID(row2, col2));
		}
	}
#undef GRID_ID(i,j)

	sort(affectedGridCell.begin(), affectedGridCell.end());
	affectedGridCell.erase(unique(affectedGridCell.begin(), affectedGridCell.end()), affectedGridCell.end());

	const int dx[8] = {-1, -1, 0, 1, 1, 1, 0, -1};
	const int dy[8] = {0, 1, 1, 1, 0, -1, -1, -1};
	//if max of sum of residual flow of 4 adj. grid cell < flow threshold, remove all edges of thegrid cell in edge choice
	for(int k = 0; k < affectedGridCell.size(); ++ k) {
		int cellId = affectedGridCell[k];
		int i = cellId / NumGridCellcol;
		int j = cellId % NumGridCellcol;

		if (i == 0 || j == 0 || i + 1 == NumGridCellrow || j + 1 == NumGridCellcol) {
			//it is on the boundary
			continue;
		}

		if (getVolcanoUpperBound(Grid, i, j) < flowThreshold) {
			Grid[i][j].noVolcano = 1;
			if(Grid[i][j].noBlackhole == 1 && Grid[i][j].noVolcano == 1) {
				for(map<int, int>::iterator it = Grid[i][j]._Edge.begin(); it != Grid[i][j]._Edge.end(); ++ it) {
					prunedCount++;
					edgeChoice[(*it).first]=1;
				}
			}
		}
	}
	fprintf(stderr, "Volcano: Dynamic Upper Bound:: pruned=%d\n",prunedCount);
}

void BlackholeinitialGridPruning(Graph &g, vector<int> &inFlow, vector<int> &outFlow, const int flowThreshold, const double areaThreshold, const int numEdge,
									vector<vector<gridCell>> &Grid, vector<int> &pickedEdge, vector<int> &edgeChoice, const vector<int> &blackholeEdge)
{
	int prunedCount = 0;
	//sum up flow of each grid
	for(int i = 0; i < Grid.size(); ++ i) {
		for(int j = 0; j < Grid[i].size(); ++ j) {
			int &sum = Grid[i][j].BlackholeresidualFlow;
			sum = 0;
			for(map<int,int>::iterator it = Grid[i][j]._Edge.begin(); it != Grid[i][j]._Edge.end(); ++ it) {
				if (blackholeEdge[it->first] != 1) {
					sum += max(inFlow[(*it).first] - outFlow[(*it).first], 0);
				}
			}
		}
	}
	//if max of sum of residual flow of 4 adj. grid cell < flow threshold, remove all edges of thegrid cell in edge choice
	for(int i = 1; i + 1 < (int)Grid.size(); ++ i) {
		for(int j = 1; j + 1 < (int)Grid[i].size(); ++ j) {
			if(getBlackholeUpperBound(Grid, i, j) < flowThreshold) {
				Grid[i][j].noBlackhole = 1;
				if(Grid[i][j].noBlackhole == 1 && Grid[i][j].noVolcano == 1) {
					for(map<int,int>::iterator it = Grid[i][j]._Edge.begin(); it != Grid[i][j]._Edge.end(); ++ it) {
						++ prunedCount;
						edgeChoice[(*it).first] = 1;
					}
				}
			}
		}
	}
	fprintf(stderr, "initialGridPruning:: pruned=%d\n", prunedCount);
}

void VolcanoinitialGridPruning(Graph& g,vector<int>& inFlow,vector<int>& outFlow,const int flowThreshold, const double areaThreshold,const int numEdge,vector<vector<gridCell>>& Grid,
						vector<int>& pickedEdge,vector<int>& edgeChoice, const vector<int> &volcanoEdge){
	int prunedCount = 0;
	//sum up flow of each grid
	for(int i = 0; i < Grid.size(); ++ i) {
		for(int j = 0; j < Grid[i].size(); ++ j) {
			int &sum = Grid[i][j].VolcanoresidualFlow;
			sum = 0;
			for (map<int,int>::iterator it = Grid[i][j]._Edge.begin(); it != Grid[i][j]._Edge.end(); ++ it) {
				if (volcanoEdge[it->first] != 1) {
					sum += max(outFlow[(*it).first] - inFlow[(*it).first], 0);
				}
			}
		}
	}
	//if max of sum of residual flow of 4 adj. grid cell < flow threshold, remove all edges of thegrid cell in edge choice
	for(int i = 1; i + 1 < (int)Grid.size(); ++ i) {
		for(int j = 1; j + 1 < (int)Grid[i].size(); ++ j) {
			if(getVolcanoUpperBound(Grid, i, j) < flowThreshold) {
				Grid[i][j].noVolcano = 1;
				if(Grid[i][j].noBlackhole == 1 && Grid[i][j].noVolcano == 1) {
					for(map<int,int>::iterator it = Grid[i][j]._Edge.begin(); it != Grid[i][j]._Edge.end(); ++ it) {
						++ prunedCount;
						edgeChoice[(*it).first] = 1;
					}
				}
			}
		}
	}
	fprintf(stderr, "initialGridPruning:: pruned=%d\n", prunedCount);
}

inline void expand(priority_queue<Entry> &candidateEdge, int u, int sign, vector<int> &mark, int stamp,
				   const vector<int> &inFlow, const vector<int> &outFlow, const vector<int> &blackholeEdge, const vector<int> &volcanoEdge,
				   set<int> &expandedNode, Graph &g, int &visitCount)
{
	if (expandedNode.count(u)) {
		return;
	}
	expandedNode.insert(u);
	vector<edge> &adjEdge = g.sg.at(u).e;
	for(int i = 0; i < adjEdge.size(); ++ i) {
		int eid = adjEdge[i].edgeID;
		int delta = inFlow[eid] - outFlow[eid];
		if(delta * sign > 0 && blackholeEdge[eid] != 1 && volcanoEdge[eid] != 1 && mark[eid] != stamp) {
			candidateEdge.push(Entry(eid, delta));
			mark[eid] = stamp;
			++ visitCount;
		}
	}
}

namespace Tarjan
{
	vector<int> low, stk;
	vector<bool> mark;
	set<pair<int, int>> vis;
	int cnt, N, n;

	void dfs(int u, const vector<vector<int>> &adj, vector<int> &color) {
		mark[u] = true;
		low[u] = ++ cnt;
		int mini = cnt;
		stk.push_back(u);
		for (int i = 0; i < adj[u].size(); ++ i) {
			int v = adj[u][i];
			if (!vis.count(make_pair(u, v))) {
				vis.insert(make_pair(u, v));
				vis.insert(make_pair(v, u));
				if (!mark[v]) {
					dfs(v, adj, color);
				}
				mini = min(mini, low[v]);
			}
		}
		if (mini == low[u]) {
			int v;
			++ N;
			do {
				v = stk.back();
				stk.pop_back();
				low[v] = n + 1;
				color[v] = N;
			} while (v != u);
		} else {
			low[u] = mini;
		}
	}

	void findBiConnect(const vector<vector<int>> &adj, vector<int> &color) {
		n = adj.size();
		vis.clear();
		mark.resize(n);
		low.resize(n);
		color.resize(n);
		stk.clear();
		for (int i = 0; i < n; ++ i) {
			mark[i] = false;
			low[i] = 0;
		}
		cnt = N = 0;
		for (int i = 0; i < n; ++ i) {
			if (!mark[i]) {
				dfs(i, adj, color);
			}
		}
	}
};

vector<int> findNonBridge(const vector<int> &edges, Graph &g)
{
	map<int, int> relabel;
	for (int i = 0; i < edges.size(); ++ i) {
		int currentEdge = edges[i];
		int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
		addNew(relabel, f);
		addNew(relabel, t);
	}
	vector<vector<int>> adj(relabel.size(), vector<int>());
	for (int i = 0; i < edges.size(); ++ i) {
		int currentEdge = edges[i];
		int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
		f = addNew(relabel, f);
		t = addNew(relabel, t);
		adj[f].push_back(t);
		adj[t].push_back(f);
	}
	for (int i = 0; i < adj.size(); ++ i) {
		sort(adj[i].begin(), adj[i].end());
		adj[i].erase(unique(adj[i].begin(), adj[i].end()), adj[i].end());
	}
	vector<int> color;
	Tarjan::findBiConnect(adj, color);

	vector<int> ret;
	for (int i = 0; i < edges.size(); ++ i) {
		int currentEdge = edges[i];
		int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
		f = addNew(relabel, f);
		t = addNew(relabel, t);
		if (color[f] == color[t]) {
			ret.push_back(currentEdge);
		}
	}
	return ret;
}

void incrementalDetection(Graph& g,vector<int>& inFlow,vector<int>& outFlow, const Configuration &config, const int numEdge,vector<vector<int>>& prevBlackholes,
						  vector<vector<int>>& prevVolcanos,vector<int>& blackholeEdge,vector<int>& volcanoEdge, int prefix_k, int prefix_i, int prefix_j, int &visitCount, vector<int> &edgeChoice)
{
	debug("incremental detection begin!");
	//for each blackhole and volcano
		//1. check the area 

		//2. if area < areathreshold, expand until area>=area threshold

		//3. if 
	vector<vector<int>> newBlackholes, newVolcanos;
	vector<int> mark(numEdge, 0);
	int stamp = 0;
	vector<vector<int>> prev = prevBlackholes;
	if (config.cross == 1) {
		for (int i = 0; i < prevVolcanos.size(); ++ i) {
			prev.push_back(prevVolcanos[i]);
		}
	}
	//Blackholes
	for (int iter = 0; iter < prev.size(); ++ iter) {
		set<int> expandedNode;
		++ stamp;
		vector<int> candidateBlackhole = prev[iter];
		priority_queue<Entry> candidateEdge; // positive & adjacent edges
		bool valid = true;
		for (int i = 0; i < candidateBlackhole.size(); ++ i) {
			int currentEdge = candidateBlackhole[i];
			if (edgeChoice[currentEdge] == 1 || blackholeEdge[currentEdge] == 1 || volcanoEdge[currentEdge] == 1) {
				valid = false;
				break;
			}
			mark[currentEdge] = stamp;
			int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
			//insert positive adj. edges of node f and t into queue
			expand(candidateEdge, f, 1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
			expand(candidateEdge, t, 1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
		}
		if (!valid) {
			continue;
		}
		const int tryTimeLimit = 10;
		for (int tryTime = 0; tryTime < tryTimeLimit && candidateEdge.size(); ++ tryTime) {
			double minLat = INT_MAX, minLong = INT_MAX, maxLat = 0, maxLong = 0;
			int absFlow = 0;
			for (int i = 0; i < candidateBlackhole.size(); ++ i) {
				int currentEdge = candidateBlackhole[i];
				absFlow += inFlow[currentEdge] - outFlow[currentEdge];
				int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
				double fLat = g.sg.at(f).la, fLong = g.sg.at(f).lo, tLat = g.sg.at(t).la, tLong = g.sg.at(t).lo;

				updateInterval(minLat, maxLat, fLat);
				updateInterval(minLat, maxLat, tLat);
				updateInterval(minLong, maxLong, fLong);
				updateInterval(minLong, maxLong, tLong);
			}
			priority_queue<Entry> nextTurn;
			while (candidateEdge.size()) {
				int currentEdge = candidateEdge.top().eid, delta = candidateEdge.top().priority;
				candidateEdge.pop();

				int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
				double fLat = g.sg.at(f).la, fLong = g.sg.at(f).lo, tLat = g.sg.at(t).la, tLong = g.sg.at(t).lo;
				if (checkAndUpdateArea(minLat, minLong, maxLat, maxLong, fLat, fLong, tLat, tLong, config.areaThreshold)) {
					candidateBlackhole.push_back(currentEdge);
					absFlow += delta;
					expand(candidateEdge, f, 1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
					expand(candidateEdge, t, 1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
				} else {
					nextTurn.push(Entry(currentEdge, delta));
				}
			}
			if (absFlow > config.flowThreshold) {
				// detected a blackhole!
				UpdateBlackholeEdge(candidateBlackhole, blackholeEdge);
				fprintf(stderr, "Incremental Blackhole Found! after %d tries\n", tryTime + 1);//blackhoel subgraph is stored in tempBlackholes
				printBlackholeToFile(candidateBlackhole, g, prefix_k, prefix_i, prefix_j, absFlow, inFlow, outFlow, config.blackholePatternFile.c_str());
				newBlackholes.push_back(candidateBlackhole);
				break;
			}
			candidateEdge = nextTurn;
			
			//delete one edge
			vector<int> canDelete = findNonBridge(candidateBlackhole, g);
			int choice = -1, mini = 0;
			for (int i = 0; i < canDelete.size(); ++ i) {
				int eid = canDelete[i], delta = inFlow[i] - outFlow[i];
				if (choice == -1 || mini > delta) {
					mini = delta;
					choice = eid;
				}
			}
			for (int i = 0; i < candidateBlackhole.size(); ++ i) {
				if (candidateBlackhole[i] == choice) {
					candidateBlackhole.erase(candidateBlackhole.begin() + i);
					//edgeChoice[choice] = 1;
					break;
				}
			}
		}
	}

	if (config.cross == 0) {
		prev = prevVolcanos;
	}
	//Volcanos
	for (int iter = 0; iter < prev.size(); ++ iter) {
		set<int> expandedNode;
		++ stamp;
		vector<int> candidateVolcano = prev[iter];
		priority_queue<Entry> candidateEdge; // positive & adjacent edges
		bool valid = true;
		for (int i = 0; i < candidateVolcano.size(); ++ i) {
			int currentEdge = candidateVolcano[i];
			if (edgeChoice[currentEdge] == 1 || blackholeEdge[currentEdge] == 1 || volcanoEdge[currentEdge] == 1) {
				valid = false;
				break;
			}
			mark[currentEdge] = stamp;
			int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
			//insert positive adj. edges of node f and t into queue
			expand(candidateEdge, f, -1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
			expand(candidateEdge, t, -1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
		}
		if (!valid) {
			continue;
		}
		const int tryTimeLimit = 10;
		for (int tryTime = 0; tryTime < tryTimeLimit && candidateEdge.size(); ++ tryTime) {
			double minLat = INT_MAX, minLong = INT_MAX, maxLat = 0, maxLong = 0;
			int absFlow = 0;
			for (int i = 0; i < candidateVolcano.size(); ++ i) {
				int currentEdge = candidateVolcano[i];
				absFlow += outFlow[currentEdge] - inFlow[currentEdge];
				int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
				double fLat = g.sg.at(f).la, fLong = g.sg.at(f).lo, tLat = g.sg.at(t).la, tLong = g.sg.at(t).lo;

				updateInterval(minLat, maxLat, fLat);
				updateInterval(minLat, maxLat, tLat);
				updateInterval(minLong, maxLong, fLong);
				updateInterval(minLong, maxLong, tLong);
			}
			priority_queue<Entry> nextTurn;
			while (candidateEdge.size()) {
				int currentEdge = candidateEdge.top().eid, delta = candidateEdge.top().priority;
				candidateEdge.pop();

				int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
				double fLat = g.sg.at(f).la, fLong = g.sg.at(f).lo, tLat = g.sg.at(t).la, tLong = g.sg.at(t).lo;
				if (checkAndUpdateArea(minLat, minLong, maxLat, maxLong, fLat, fLong, tLat, tLong, config.areaThreshold)) {
					candidateVolcano.push_back(currentEdge);
					absFlow += delta;
					expand(candidateEdge, f, -1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
					expand(candidateEdge, t, -1, mark, stamp, inFlow, outFlow, blackholeEdge, volcanoEdge, expandedNode, g, visitCount);
				} else {
					nextTurn.push(Entry(currentEdge, delta));
				}
			}
			if (absFlow > config.flowThreshold) {
				// detected a blackhole!
				UpdateBlackholeEdge(candidateVolcano, volcanoEdge);
				fprintf(stderr, "Incremental Volcano Found! after %d tries\n", tryTime + 1);
				printVolcanoToFile(candidateVolcano, g, prefix_k, prefix_i, prefix_j, absFlow, inFlow, outFlow, config.volcanoPatternFile.c_str());
				newVolcanos.push_back(candidateVolcano);
				break;
			}
			candidateEdge = nextTurn;
			
			//delete one edge
			vector<int> canDelete = findNonBridge(candidateVolcano, g);
			int choice = -1, mini = 0;
			for (int i = 0; i < canDelete.size(); ++ i) {
				int eid = canDelete[i], delta = outFlow[i] - inFlow[i];
				if (choice == -1 || mini > delta) {
					mini = delta;
					choice = eid;
				}
			}
			for (int i = 0; i < candidateVolcano.size(); ++ i) {
				if (candidateVolcano[i] == choice) {
					candidateVolcano.erase(candidateVolcano.begin() + i);
					//edgeChoice[choice] = 1;
					break;
				}
			}
		}
	}

	prevBlackholes = newBlackholes;
	prevVolcanos = newVolcanos;

	debug("incremental detection finishied!");
}

const int maxlength = 10000000;
char line[maxlength + 1];

inline void tryLoad(vector<vector<int>> &prev, char* filename)
{
	FILE* f = fopen(filename, "r");
	if (f != NULL) {
		while (true) {
			for (;fgets(line, maxlength, f) && line[0] != 'c';);
			if (line[0] != 'c') {
				break;
			}
			vector<int> subgraph;
			for (;fgets(line, maxlength, f) && isdigit(line[0]);) {
				stringstream in(line);
				int id ,edge;
				in >> id >> edge;
				subgraph.push_back(edge);
			}
			prev.push_back(subgraph);
		}
		fclose(f);
	} else {
		fprintf(stderr, "%s not found!\n", filename);
	}
}

void historicalDetection(vector<vector<int>> &prevB, vector<vector<int>> prevV, const Configuration &config, int prefix_k, int day, int period) {
	//... not yet implemented
	// load historical blackholes and volcanos
	for (int preDay = day - 7; preDay > 0; preDay -= 7) {
		char filename[255];
		sprintf(filename, "%s%d_%d-%d.blackhold", config.historicalFile.c_str(), preDay, period, period + config.stride - 1);
		tryLoad(prevB, filename);

		sprintf(filename, "%s%d_%d-%d.volcano", config.historicalFile.c_str(), preDay, period, period + config.stride - 1);
		tryLoad(prevV, filename);
	}
	fprintf(stderr, "%d + %d loaded!\n", prevB.size(), prevV.size());
}


void continuousDetection(Graph &g, vector<int> &inFlow, vector<int> &outFlow,const Configuration &config, const int numEdge,
						 vector<vector<int>> &prevBlackholes, vector<vector<int>> &prevVolcanos, vector<int> &blackholeEdge,vector<int> &volcanoEdge,
						 int prefix_k, int prefix_i, int prefix_j, int &visitCount, vector<int> &edgeChoice)
{
	if (config.noPrev == 1) {
		prevBlackholes.clear();
		prevVolcanos.clear();
	}
	if (config.noHis == 0) {
		historicalDetection(prevBlackholes, prevVolcanos, config, prefix_k, prefix_i, prefix_j);
	}
	incrementalDetection(g,inFlow,outFlow,config,numEdge,prevBlackholes,prevVolcanos,blackholeEdge,volcanoEdge, prefix_k, prefix_i, prefix_j, visitCount, edgeChoice);
}

void RemoveCurrentSubgraphfromEdgeChoice(vector<int> &currentSubgraph, vector<int> &edgeChoice)
{
	for(int i = 0; i < currentSubgraph.size(); ++ i) {
		edgeChoice[currentSubgraph[i]] = 1;
	}
}

void sortEdgeFlowIndecsendingOrder(vector<data> &Num2, int numEdge, vector<int> &inFlow, vector<int> &outFlow, const char* absFlowFile)
{
	FILE * pFile = fopen(absFlowFile, "w");
	for(int i = 0; i < numEdge; ++ i) {
		data d;
		d.edgeID = i;
		d.absFlow = inFlow[i] - outFlow[i];
		Num2.push_back(d);
		fprintf(pFile,"%d %d\n", d.edgeID, d.absFlow);
	}
	fclose(pFile);
	std::sort(Num2.begin(), Num2.end(), compareByDescAbsFlow);
}

void ShuffleEdgeRandom(vector<int> &order, int numEdge)
{
	order.resize(numEdge);
	for(int i = 0; i < numEdge; ++ i) {
		order[i] = i;
	}
	random_shuffle(order.begin(), order.end());
}

bool tryPushEdge(int eid, vector<int> &pickedEdge, vector<int> &blackholeEdge, vector<int> &volcanoEdge, vector<int> &edgeChoice, int &SkipCount)
{
	if(pickedEdge[eid] == 1 || blackholeEdge[eid] == 1 || volcanoEdge[eid] == 1 || edgeChoice[eid]==1) {
		++ SkipCount;
		return true;
	}
	return false;
}

void cleanUp(int numEdge, int &absInFlow, int &absOutFlow, const Configuration &config,
			 int &blackholeCount, int &tempBlackholesFlow, int &tempVolcanosFlow, int &volcanoCount,
			 int &visitCount, int &totalVisitCount,
			 int &forloopCount, double &area,
			 vector<int> &tempBlackholes, vector<int> &tempVolcanos, vector<int> &blackholeEdge, vector<int> &edgeChoice,
			 double &maxLat, double &maxLong, double &minLat, double &minLong,
			 vector<vector<int>> &prevBlackholes, vector<vector<int>> &prevVolcanos,
			 int useDynamicUpperBound, vector<int> &inFlow, vector<int> &outFlow,
			 vector<vector<gridCell>> &Grid, vector<int> &currentSubgraph,
			 vector<int> &pickedEdge, vector<int> &volcanoEdge, vector<int> &visited, vector<int> &visitedEdge,
			 Graph &g, int k, int i, int j)
{
	//clean up after expanding from an initial edge
	if(absInFlow > config.flowThreshold) {
		// a blackhole detected!
		UpdateBlackholeEdge(tempBlackholes, blackholeEdge);
		RemoveBlackholeEdgefromEdgeChoice2(tempBlackholes, edgeChoice);
//		fprintf(stderr, "Blackhole Found!\n");//blackhoel subgraph is stored in tempBlackholes
//		fprintf(stderr, "    edge visited %d currentSubgraph total flow %d\n", visitCount, tempBlackholesFlow);
		++ blackholeCount;

		printBlackholeToFile(tempBlackholes, g, k, i, j, tempBlackholesFlow, inFlow, outFlow, config.blackholePatternFile.c_str());

		//dynamic upper bound technique for blackhole
		if(useDynamicUpperBound==1) {
			blackholedynamicallyComputeUpperBound(g, inFlow, outFlow, config.flowThreshold, config.areaThreshold, numEdge, Grid, pickedEdge, edgeChoice, tempBlackholes, maxLat, maxLong, minLat, minLong);
		}
		prevBlackholes.push_back(tempBlackholes);
	} else if (absOutFlow > config.flowThreshold) {
		// a volcano detected!
		UpdateBlackholeEdge(tempVolcanos, volcanoEdge);
		RemoveBlackholeEdgefromEdgeChoice2(tempVolcanos, edgeChoice);
//		fprintf(stderr, "Volcano Found!\n");//blackhoel subgraph is stored in tempBlackholes
//		fprintf(stderr, "    edge visited %d currentSubgraph total flow %d\n",visitCount,tempVolcanosFlow);
		++ volcanoCount;

		printVolcanoToFile(tempVolcanos, g, k, i, j, tempVolcanosFlow, inFlow, outFlow, config.volcanoPatternFile.c_str());

		//dynamic upper bound technique for volcano
		if(useDynamicUpperBound==1) {
			volcanodynamicallyComputeUpperBound(g, inFlow, outFlow, config.flowThreshold, config.areaThreshold, numEdge, Grid, pickedEdge, edgeChoice, tempVolcanos, maxLat, maxLong, minLat, minLong);
		}
		prevVolcanos.push_back(tempVolcanos);
	} else {
		//fprintf(stderr, "No Blackhole or Volcano Found!\n");
		//fprintf(stderr, "edge visited %d currentSubgraph total inflow %d outflow %d\n",visitCount,absInFlow,absOutFlow);
	}
	currentSubgraph.clear();
	tempBlackholes.clear();
	tempVolcanos.clear();
//	fprintf(stderr, "forloopCount = %d\n", forloopCount);

	totalVisitCount += visitCount;
	tempBlackholesFlow = tempVolcanosFlow = absInFlow = absOutFlow = visitCount = area = forloopCount = 0;
	for(int i=0;i<visited.size();i++) {
		visitedEdge[visited[i]]=0;
	}
	visited.clear();
//	fprintf(stderr, "\n");
}

inline double getPriority(int iter, int eid, const vector<int> &inFlow, const vector<int> &outFlow, int useExpansionOrdering)
{
	if (useExpansionOrdering) {
		if (iter == 0) {
			return + inFlow[eid] - outFlow[eid];
		} else {
			return - inFlow[eid] + outFlow[eid];
		}
	} else {
		return 0;
	}
}

//this is the vanilla baseline
int VolcanoBlackholeDetection(Graph &g, vector<int> &inFlow, vector<int> &outFlow, const Configuration &config,
							  const int numEdge, vector<vector<gridCell>> &Grid,
							  double maxLat, double maxLong, double minLat, double minLong,
							  vector<vector<int>> &prevBlackholes, vector<vector<int>>& prevVolcanos,
							  int k2, int i, int j, int &continuousBlackhole, int &continuousVolcano)
{
	if (config.dump == 1) {
		vector<int> pos, neg;
		int ps = 0, ns = 0;
		for (int i = 0; i < numEdge; ++ i) {
			int delta = inFlow[i] - outFlow[i];
			if (delta > 0) {
				pos.push_back(i);
				ps += delta;
			} else if (delta < 0) {
				neg.push_back(i);
				ns -= delta;
			}
		}
		printBlackholeToFile(pos, g, k2, i, j, ps, inFlow, outFlow, config.blackholePatternFile.c_str());
		prevBlackholes.clear(); prevBlackholes.push_back(pos);

		printVolcanoToFile(neg, g, k2, i, j, ns, inFlow, outFlow, config.volcanoPatternFile.c_str());
		prevVolcanos.clear(); prevVolcanos.push_back(neg);
		return 0;
	}
	//declare parameters
	vector<int> pickedEdge,//store edges that have been used as initial edge
		visitedEdge,//store edges that have been visited in 1 exoansion, reset every iteration
		blackholeEdge,//store edges that are blackhole
		volcanoEdge,//store edges that are volcano
		visited;//this one store the visited edge for each expansion. for restoring the visitedEdge array (this one no need to initialize!)
	priority_queue<Entry> q; //a pirority queue sort based on our function
	vector<int> edgeChoice, currentSubgraph, order;
	vector<int> tempBlackholes, tempVolcanos;
	double area = 0;
	int blackholeCount=0,volcanoCount=0,totalVisitCount=0,visitCount=0,forloopCount=0,absInFlow=0,absOutFlow=0,tempBlackholesFlow=0,tempVolcanosFlow=0,SkipCount=0;

	FILE * pFile2=fopen(config.infoFile.c_str(),"a");
	int startTime, endTime, totalTime=0;//*** timer
	startTime = time(NULL);//*** timer

	//initialization
	pickedEdge.assign(numEdge,0); visitedEdge.assign(numEdge,0); blackholeEdge.assign(numEdge,0); volcanoEdge.assign(numEdge,0); edgeChoice.assign(numEdge,0);

	//continuous blackhole detection //not yet implemented!
	if(config.useContinuous == 1) {
		continuousDetection(g, inFlow, outFlow, config, numEdge,prevBlackholes,prevVolcanos,blackholeEdge,volcanoEdge,k2,i,j, visitCount, edgeChoice);
		blackholeCount = prevBlackholes.size();
		volcanoCount = prevVolcanos.size();
	} else {
		prevBlackholes.clear();
		prevVolcanos.clear();
	}

	continuousBlackhole = prevBlackholes.size();
	continuousVolcano = prevVolcanos.size();

	//initial grid cell pruning //not yet implemented!
	if(config.useInitialPruning == 1 || config.useDynamicUpperBound==1) {
		BlackholeinitialGridPruning(g,inFlow,outFlow,config.flowThreshold,config.areaThreshold,numEdge,Grid,pickedEdge,edgeChoice, blackholeEdge);
		VolcanoinitialGridPruning(g,inFlow,outFlow,config.flowThreshold,config.areaThreshold,numEdge,Grid,pickedEdge,edgeChoice, volcanoEdge);
	}

	if(config.useInitEdgePickOrdering == 1) {
		vector<data> temp;
		sortEdgeFlowIndecsendingOrder(temp,numEdge,inFlow,outFlow, config.absFlowFile.c_str());
		order.resize(numEdge);
		for (int i = 0; i < numEdge; ++ i) {
			order[i] = temp[i].edgeID;
		}
	} else {
		ShuffleEdgeRandom(order,numEdge);
	}

	for(int k = 0; k < numEdge; ++ k) {
//////////////////////jsut for testing only////////////////////////////////////////////
		if(k % 40000 == 0) {
			FILE * pFile3=fopen(config.timeFile.c_str(),"a");
			endTime = time(NULL);//*** timer
			fprintf(pFile3,"edge %d %d\n", k, endTime - startTime);
			startTime = time(NULL);//*** timer
			fclose(pFile3);
		}
//////////////////////////////////////////////////////////////////////////////////////
		for (int iter = 0; iter < 2; ++ iter) {
			//one for blackhole, the other for volcano
			if (true) {
				int eid = order[k];
				if(tryPushEdge(eid,pickedEdge,blackholeEdge,volcanoEdge,edgeChoice,SkipCount)) {
					//this edge is used before
					break;
				} else {
					q.push(Entry(eid, getPriority(iter, eid, inFlow, outFlow, config.useExpansionOrdering)));
					visited.push_back(eid);
					visitedEdge[eid] = 1;
					++ visitCount;
					//fprintf(stderr, "bfs time = %d, initialEdgeChoice %d\n", k, eid);
				}
			}

			double BH_minLat = INT_MAX, BH_minLong = INT_MAX, BH_maxLat = 0, BH_maxLong = 0;
			bool island = true;
			//BFS (use queue) expansion from initial edge until area==area threshold, check blackhole
			while (q.size()) {
				//insert an edge into queue if the edge connect to either node of the poped edge
				int currentEdge = q.top().eid;
				q.pop();

				int delta = inFlow[currentEdge]-outFlow[currentEdge];

				if (config.useExpansionOrdering == 1) {
					if (iter == 0 && delta < 0 || iter == 1 && delta > 0) {
						while (q.size()) {
							q.pop();
						}
						break;
					}
				}


				//retrieve the f node and t node of the currentEdge
				int f = g.eg.at(currentEdge).nodef, t = g.eg.at(currentEdge).nodet;
				double fLat = g.sg.at(f).la, fLong = g.sg.at(f).lo, tLat = g.sg.at(t).la, tLong = g.sg.at(t).lo;

				//check area; area is updated if return true; if return false, area is not updated
				bool within = checkAndUpdateArea(BH_minLat, BH_minLong, BH_maxLat, BH_maxLong, fLat, fLong, tLat, tLong, config.areaThreshold);

				if (within) {
					//update currentSubraph
					currentSubgraph.push_back(currentEdge);
					//compute flow
					absInFlow += delta;
					absOutFlow += -delta;

					//update tempBlackholes if discover a blackhole
					if (absInFlow > config.flowThreshold) {
						tempBlackholes = currentSubgraph;
						tempBlackholesFlow = absInFlow;
					}
					//update tempVolcanos if discover a volcano
					if (absOutFlow > config.flowThreshold) {
						tempVolcanos = currentSubgraph;
						tempVolcanosFlow = absOutFlow;
					}
					//insert adj. edges of node f and t into queue
					if (true) {
						// f side
						vector<edge> &adjEdge = g.sg.at(f).e;
						for(int i = 0; i < adjEdge.size(); ++ i) {
							forloopCount ++;
							int eid = adjEdge[i].edgeID;
							++ visitCount;
							if(blackholeEdge[eid] != 1 && volcanoEdge[eid] != 1 && visitedEdge[eid] != 1) {
								q.push(Entry(eid, getPriority(iter, eid, inFlow, outFlow, config.useExpansionOrdering)));
								visitedEdge[eid] = 1;
								visited.push_back(eid);
							}
						}
					}
					if (true) {
						// t side
						vector<edge> &adjEdge = g.sg.at(t).e;
						for(int i = 0; i < adjEdge.size(); ++ i) {
							forloopCount ++;
							int eid = adjEdge[i].edgeID;
							++ visitCount;
							if(blackholeEdge[eid] != 1 && volcanoEdge[eid] != 1 && visitedEdge[eid] != 1) {
								q.push(Entry(eid, getPriority(iter, eid, inFlow, outFlow, config.useExpansionOrdering)));
								visitedEdge[eid] = 1;
								visited.push_back(eid);
							}
						}
					}
				} else {
					island = false;
				}
			}
		
			if(island == true && config.useIslandRemoval == 1) {
				RemoveCurrentSubgraphfromEdgeChoice(currentSubgraph, edgeChoice);
			}

			cleanUp(numEdge, absInFlow, absOutFlow, config, blackholeCount, tempBlackholesFlow, tempVolcanosFlow, volcanoCount, visitCount, totalVisitCount, forloopCount, area
				 ,tempBlackholes,tempVolcanos,blackholeEdge,edgeChoice,maxLat,maxLong,minLat,minLong
				 ,prevBlackholes,prevVolcanos,config.useDynamicUpperBound,inFlow,outFlow
				 ,Grid,currentSubgraph,pickedEdge,volcanoEdge,visited,visitedEdge,g,k2,i,j);
		}
	}
	fprintf(pFile2,"SkipCount= %d blackholeCount=%d VolcanoCount=%d\n",SkipCount,blackholeCount,volcanoCount);
	fclose(pFile2);
	fprintf(stderr, "No. of blackhole Found %d\n totalVisitCount=%d\n",blackholeCount,totalVisitCount);
	fprintf(stderr, "No. of volcano Found %d\n \n",volcanoCount);

	return totalVisitCount;
}


#endif