#include <vector>
#include <algorithm>
#include <utility>
#include <numeric>
#include <map>
#include <iostream>
#include <cmath>
#include <map>
// #include <sstream>  // if c++ 11 is not available
#include <string>

// #include "boost/tuple/tuple.hpp"
// #include <tuple>

#include "TH2D.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TString.h"
#include "TGraph.h"
#include "TEllipse.h"
#include "TF1.h"
#include "TMath.h"

#include "PndAhTools.h"
#include "PndAhHoughTransform.h"
// #include "PndAhMultipeakAnalysis.h"
#include "PndAhConformalMap.h"

class PndAhIterativeMaximumDeleter {
public:
	PndAhIterativeMaximumDeleter();
	PndAhIterativeMaximumDeleter(houghInfo & info, unsigned int iterations = 1) : fVX(info.x), fVY(info.y), fVR(info.r), fMaxAngle(info.maxAngle), fEveryXDegrees(info.everyXDegrees), fVerbose(info.verbose), fIterations(iterations) {DoEverthing();};
	PndAhIterativeMaximumDeleter(std::vector<double> x, std::vector<double> y, std::vector<double> r, double maxAngle = 360, double everyXDegrees = 30, bool verbose = false) : fVX(x), fVY(y), fVR(r), fMaxAngle(maxAngle), fEveryXDegrees(everyXDegrees), fVerbose(verbose) {DoEverthing();};
	~PndAhIterativeMaximumDeleter();

	void SetIterations(unsigned int iterations = 1) {fIterations = iterations;};

	std::vector<int> GetMaxBins() {return fVMaxBins;};
	std::vector<std::pair<double, double> > GetMaxCoordinates() {return fVMaxCoordinates;}
	std::vector<TH2D *> GetHoughHistograms() {return fVHoughHists;};
	std::vector<std::pair<double, double> > GetBorderXCoordinates() {return fBorderXCoordinates;};
protected:
	void DoEverthing() {
		PrintEverything();
		std::cout << " -- Before Iterating -- " << std::endl;
		for (int i = 0; i < fIterations; ++i) {
			std::cout << " -- Starting Iteration " << i << std::endl;
			if (PndAhHoughTransform::vectorsAreEmpty(fVX, fVY, fVR)) {
				break;
			} else {
				PndAhHoughTransform * hough = new PndAhHoughTransform(fVX, fVY, fVR, fMaxAngle, fEveryXDegrees, fVerbose, true);
				// std::cout << " -- Getting Histogram." << std::endl;
				TH2D * houghHist = hough->GetHistogram(true, "houghIt" + PndAhTools::numberToString(i));

				if (stopIterating(houghHist)) break;  // some stop criterion is reached
				// std::cout << " -- Getting Max Passings." << std::endl;
				std::vector<unsigned int> passingHits = generateListOfMaxPassingHits(hough, houghHist);

				// std::cout << " -- passingHits.size() = " << passingHits.size() << std::endl;

				sortMaxPassingHits(passingHits);

				removeMaxPassingHits(passingHits);

				fVHoughHists.push_back(houghHist);
				fVMaxBins.push_back(getHistMaxBinNumber(houghHist));
				fVMaxCoordinates.push_back(getHistMaximumCoordinates(houghHist));
				// std::cout << " -- Ending Iteration " << i << std::endl;
			}
		}
		// std::cout << " -- After Iterating -- " << std::endl;
	}
	bool stopIterating(TH2D * hist) {
		if (stopCriterion_histIsBelowZThreshold(hist, 3)) return true;
		return false;
	}
	bool stopCriterion_histIsBelowZThreshold(TH2D * hist, double threshold) {
		return hist->GetMaximum() <= threshold;
	}
	void PrintEverything() {
		std::cout << "fVX.size() = " << fVX.size() << std::endl;
		std::cout << "fMaxAngle = " << fMaxAngle << std::endl;
		std::cout << "fEveryXDegrees = " << fEveryXDegrees << std::endl;
		std::cout << "fVerbose = " << fVerbose << std::endl;
		if (fIterations) std::cout << "fIterations = " << fIterations << std::endl;
	}
	int getHistMaxBinNumber(TH2D * hist) {
		return hist->GetMaximumBin();
	}
	std::pair<double, double> getHistMaximumCoordinates(TH2D * hist) {
		int x, y, z;
		hist->GetMaximumBin(x, y, z);

		double xCoordinate = hist->GetXaxis()->GetBinCenter(x);
		double yCoordinate = hist->GetYaxis()->GetBinCenter(y);

		return std::make_pair(xCoordinate, yCoordinate);
	}
	std::vector<unsigned int> generateListOfMaxPassingHits(PndAhHoughTransform * ht, TH2D * houghHist) {
		std::multimap<int, unsigned int> binIdsAndHitIds = ht->GetBinIdHitId();
		int maxBinId = getHistMaxBinNumber(houghHist);

		// std::cout << " -- maxBinId = " << maxBinId << std::endl;
			// if (houghHist->GetBinContent(maxBinId) < 0.0001) {  // this only works as long as there are no negative bin contents!
			// 	std::cout << " -- It's the end of the world!" << std::endl;
			// 	fHistIsEmpty = true;
			// 	std::vector<unsigned int> emptyVector;
			// 	return emptyVector;
			// }

		std::multimap<int, unsigned int>::iterator itLow, itUp;
		itLow = binIdsAndHitIds.lower_bound(maxBinId);
		itUp = binIdsAndHitIds.upper_bound(maxBinId);

		std::vector<unsigned int> vOfGoodIds;

		for (std::multimap<int, unsigned int>::iterator it = itLow; it != itUp; ++it) {
			vOfGoodIds.push_back((*it).second);
		}
		return vOfGoodIds;
	}
	void sortMaxPassingHits(std::vector<unsigned int> & list) {
		// PndAhTools::printVector(list);
		std::sort(list.begin(), list.end());
		// PndAhTools::printVector(list);
		std::reverse(list.begin(), list.end());
		// PndAhTools::printVector(list);
	}
	void removeMaxPassingHits(std::vector<unsigned int> list) {
		double tempXInner, tempXOuter;
		bool firstLoop = true;
		for (std::vector<unsigned int>::iterator it = list.begin(); it < list.end(); it++) {
			double currentXvalue = fVX.at((*it));
			// std::cout << " -- removeMaxPassingHits -- currentXvalue = " << currentXvalue << std::endl;
			if (firstLoop) {
				// std::cout << " -- removeMaxPassingHits -- firstLoop" << std::endl;
				tempXInner = currentXvalue;
				tempXOuter = currentXvalue;
				firstLoop = false;
			}
			// std::cout << " -- removeMaxPassingHits -- abs(tempXInner), abs(tempXOuter), abs(currentXvalue)" << std::abs(tempXInner) << ", " << std::abs(tempXOuter) << ", " << std::abs(currentXvalue) << std::endl;
			if (std::abs(currentXvalue) < std::abs(tempXInner)) tempXInner = currentXvalue;
			if (std::abs(currentXvalue) > std::abs(tempXOuter))	tempXOuter = currentXvalue;

			// std::cout << " -- (*it) = " << (*it) << std::endl;
			// std::cout << " -- fVX.size() = " << fVX.size() << std::endl;
			fVX.erase(fVX.begin() + (*it));
			fVY.erase(fVY.begin() + (*it));
			fVR.erase(fVR.begin() + (*it));
			// std::cout << " -- removed for (*it) = " << (*it) << std::endl;
		}
		// std::cout << " -- removeMaxPassingHits -- " << tempXInner << ", " << tempXOuter << std::endl;
		fBorderXCoordinates.push_back(std::make_pair(tempXInner, tempXOuter));
		// std::cout << " -- End of removeMaxPassingHits()" << std::endl;
	}
private:
	std::vector<double> fVX;
	std::vector<double> fVY;
	std::vector<double> fVR;
	double fMaxAngle;
	double fEveryXDegrees;
	double fVerbose;

	unsigned int fIterations;

	std::vector<int> fVMaxBins;
	std::vector<std::pair<double, double> > fVMaxCoordinates;
	std::vector<TH2D *> fVHoughHists;
	std::vector<std::pair<double, double> > fBorderXCoordinates;  // inner, outer

	// bool fHistIsEmpty;
	
};

TGraph vectorsToGraph(std::vector<double> x, std::vector<double> y) {
	TGraph vectorGraph;
	for (int i = 0; i < x.size(); ++i) {
		vectorGraph.SetPoint(i, x[i], y[i]);
	}
	vectorGraph.SetMarkerStyle(31);
	return vectorGraph;
}

TF1 * houghParametersToLine(double r, double alpha) {
	TF1 * equation = new TF1("houghLine", "[0]/sin([1]) - x*cos([1])/sin([1])");
	equation->SetParameter(0, r);
	equation->SetParameter(1, alpha * TMath::DegToRad());
	return equation;
}

std::vector<TEllipse *> vectorsToEllipses(std::vector<double> x, std::vector<double> y, std::vector<double> r) {
	std::vector<TEllipse *> ellipses;
	for (int i = 0; i < x.size(); ++i) {
		TEllipse * thisEllipse = new TEllipse(x[i], y[i], r[i], r[i]);
		thisEllipse->SetFillStyle(0);
		ellipses.push_back(thisEllipse);
	}
	return ellipses;
}

std::vector<TF1 *> maximaToTracks(std::vector<std::pair<double, double> > vOfMaxima, std::vector<std::pair<double, double> > vOfBorders) {
	std::vector<TF1 *> tempVector;

	// std::vector<std::pair<double, double> >::iterator itMax = vOfMaxima.begin();
	// std::vector<std::pair<double, double> >::iterator itBorder = vOfBorders.begin();
	for (std::vector<std::pair<double, double> >::iterator itMax = vOfMaxima.begin(), itBorder = vOfBorders.begin(); itMax != vOfMaxima.end() && itBorder != vOfBorders.end(); itMax++, itBorder++) {  // initializing two distinct iterators here only works because they are of the same type - if they are not, the iterators have to be initialized outside the for loop individually. another way of taking care of this kind of thing is using zip operators
		double currentMaxAlpha = std::get<0>(*itMax);
		double currentMaxRho = std::get<1>(*itMax);
		TF1 * currentTrack = houghParametersToLine(currentMaxRho, currentMaxAlpha);

		double currentInnerLimit = std::get<0>(*itBorder) * 0.95;
		double currentOuterLimit = std::get<1>(*itBorder) * 1.05;
		currentTrack->SetRange(currentInnerLimit, currentOuterLimit);

		tempVector.push_back(currentTrack);
	}
	return tempVector;
}


int main(int argc, char** argv) {
	double everyXDegrees = 1;
	if (argc > 1) everyXDegrees = static_cast<double>(atof(argv[1]));
	int readInUpToLineNumber = 20;
	if (argc > 2) readInUpToLineNumber = static_cast<int>(atof(argv[2]));

	// std::vector<double> origX, origY, origR;

	// PndAhTools::readPoints("real_data.txt", origX , origY, origR, readInUpToLineNumber);

	// PndAhConformalMap * confMap = new PndAhConformalMap(origX, origY, origR);

	// std::vector<double> x = confMap->GetConfMappedXVector();
	// std::vector<double> y = confMap->GetConfMappedYVector();
	// std::vector<double> r = confMap->GetConfMappedRVector();

	std::vector<double> x, y, r;

	PndAhTools::readPoints("real_data.txt", x , y, r, readInUpToLineNumber);

	std::vector<double> fakeR(x.size());
	std::fill(fakeR.begin(), fakeR.end(), 0.);

	unsigned int maximumIterations = 8;

	houghInfo currentHoughInfo;
	currentHoughInfo.x = x;
	currentHoughInfo.y = y;
	currentHoughInfo.r = fakeR;
	currentHoughInfo.maxAngle = 180.;
	currentHoughInfo.everyXDegrees = everyXDegrees;
	currentHoughInfo.verbose = false;

	PndAhIterativeMaximumDeleter * myDeleter = new PndAhIterativeMaximumDeleter(currentHoughInfo, maximumIterations);

	auto maxBins = myDeleter->GetMaxBins();
	std::vector<TH2D *> houghHists = myDeleter->GetHoughHistograms();

	unsigned int nHist = houghHists.size();

	TApplication * theApp = new TApplication("app", &argc, argv, 0, -1);
	TCanvas * cHough = new TCanvas("cHough","cHough",5,5,1000,800);
	PndAhTools::ROOT::divideCanvas(cHough, (int)nHist, 2);
	for (int i = 0; i < nHist; ++i) {
		cHough->cd(i + 1);
		houghHists[i]->Draw("COLZ");
	}

	TCanvas * cPoints = new TCanvas("cPoints", "cPoints", 0, 0, 500, 500);
	TGraph hitPoints = vectorsToGraph(x, y);
	hitPoints.Draw("AP");
	std::vector<TEllipse *> hitEllipses = vectorsToEllipses(x, y, r);
	for(std::vector<TEllipse *>::iterator it = hitEllipses.begin(); it != hitEllipses.end(); it++) {
		(*it)->Draw("SAME");
	}

	auto maxCoordinates = myDeleter->GetMaxCoordinates();
	auto borderXCoordinates = myDeleter->GetBorderXCoordinates();
	std::vector<TF1 *> vOfTracks = maximaToTracks(maxCoordinates, borderXCoordinates);

	Color_t someColors[] = {kBlue, kGreen, kRed, kOrange, kMagenta};
	int tempI = 0;
	for(std::vector<TF1 *>::iterator it = vOfTracks.begin(); it != vOfTracks.end(); it++) {
		(*it)->SetLineColor(someColors[tempI++]);
		(*it)->Draw("SAME");
	}

	theApp->Run();
	return 1;
}