text:
/*
 * class: signalFirstInternal.cpp
 *
 * sgeisselsoeder
 *
 * Date Okt 2013
 *
 * (c) 2013 Antares Collaboration
 */

#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <numeric>
#include <omp.h>
#include <stdlib.h>
#include <string>
#include <time.h>
#include <vector>
#include <boost/lexical_cast.hpp>

#include "signalFirstInternal.h"
#include "allScales.h"
#include "candidate.hpp"
#include "clusters.h"
#include "rememberBoostHistogram.h"
#include "utils.hpp"
#include "vtkWriter.h"

#include "sgcu/sgcu.hpp"

namespace SG{

void signalFirstInternal::setOutputFolder(std::string newOutputFolder)
{
	outputFolder_ = newOutputFolder;
}

void addSourceGaussian(std::vector<SG::candidate_t>& artificial, double az, double zen, double sigma, unsigned int num){
	unsigned int counter = 0;
	double sigMax = 3.0*sigma;
	candidate_t center(az, zen);
	while (counter < num)
	{
		double offAz = sgcu::goodRand(2.0*sigMax)-sigMax;
		double offZen = sgcu::goodRand(2.0*sigMax)-sigMax;
		candidate_t temp(az+offAz, zen+offZen);

		double distance = distanceRad(center, temp);
		double sig = distance/sigma;
		double prob = std::exp(-0.5*sig*sig);

		double keepIt = sgcu::goodRand(1.0);
		if (keepIt <= prob)
		{
			artificial.push_back(temp);
			++counter;
		}
	}
}

void addSourceTemplate(std::vector<SG::candidate_t>& artificial, SG::sphericalGrid<double> const & sourceShape, unsigned int num)
{
	unsigned int counter = 0;
	unsigned int counter2 = 0;
	// TODO: check that the shape isn't empty

	// TODO: normalize template to 0:1

	while (counter < num)
	{
		++counter2;
		double offAz = sgcu::goodRand(2.0*M_PI)-M_PI;
		double offZen = sgcu::goodRand(2.0*M_PI)-M_PI;
		candidate_t temp(offAz, offZen);

		// TODO: keep the candidate according to the probability at its closest pixel in the source template
		double prob = sourceShape.findGridpoint(offAz, offZen)->proximity_;
		double keepIt = sgcu::goodRand(1.0);
		if (keepIt <= prob)
		{
			artificial.push_back(temp);
			++counter;
		}
		if (counter2 > 10000*num)
		{
			std::cerr << "DEBUG: addSourceTemplate needs insanely many tries to fill the source template. Something seems to be wrong? Aborting the computation here." << std::endl;
			throw ("party");
		}
	}
}

void addSource(std::vector<candidate_t>& artificial, double az, double zen, double deltaAz, double deltaZen, unsigned int num){
	for (unsigned int i = 0; i < num; ++i){
		double d = sgcu::goodRand(2.0*deltaAz)-deltaAz;
		double d2 = sgcu::goodRand(2.0*deltaZen)-deltaZen;
		artificial.push_back(candidate_t(az+d, zen+d2));
	}
}

void addSource(std::vector<candidate_t>& artificial, double az, double zen, double delta, unsigned int num)
{
	addSource(artificial, az, zen, delta, delta, num);
}

void dumpEvents(std::vector<candidate_t> const & events, std::string outputFile){
	std::cout << "dumpEvents: dumping " << events.size() << " events to file " << outputFile << std::endl;
	std::ofstream raus(outputFile.c_str());
	for (unsigned int i = 0; i < events.size(); ++i)
	{
		raus << events[i] << std::endl;
	}
	raus.close();
}

void dumpDebugOutputOnZeroOrTwenty(unsigned int n, std::vector<double> data, std::string const & filenameWithoutNumberOrEnding)
{
	if ( (n==0) || (n==20) )
	{
//		sgcu::saveVectorToFile(data, std::string(filenameWithoutNumberOrEnding + std::to_string(n) + ".txt") );
		sgcu::saveVectorToFile(data, std::string(filenameWithoutNumberOrEnding + boost::lexical_cast<std::string>(n) + ".txt") );
	}
}

// TODO substitute this?
void signalFirstInternal::writeNormalizationGrid(unsigned int distanceBinToPlot, std::string const & filenameWithoutNumberOrEnding)
{
	for (unsigned int i = 0; i < ratios_.aroundRatioAll.size(); ++i)
		grid_[i].proximity_ = ratios_.aroundRatioAll[i][distanceBinToPlot];
//	myVtkWriter_.write(grid_.gridpoints_, std::string(filenameWithoutNumberOrEnding + std::to_string(distanceBinToPlot) + ".vtk"));
	myVtkWriter_.write(grid_.gridpoints_, std::string(filenameWithoutNumberOrEnding + boost::lexical_cast<std::string>(distanceBinToPlot) + ".vtk"));
}

void signalFirstInternal::writeNormalizationGrids(std::string folderPrefix)
{
	writeNormalizationGrid(0, folderPrefix+"gridNormalizationNew");
	writeNormalizationGrid(20, folderPrefix+"gridNormalizationNew");
	writeNormalizationGrid(100, folderPrefix+"gridNormalizationNew");
}

void computeRadialSymmetricValuesOneScale(allScales& ratios, SG::sphericalGrid<double>& grid, unsigned int scale, std::vector<double>& radialSymmetricValues, double& allSum)
{
	//double allSum = 0.0;
	//std::vector<double> radialSymmetricValues;
	std::vector<unsigned int> counters;
	double currentDec = -99999999999.9;
	// Now store the normalization ratios on a grid
	for (unsigned int i = 0; i < ratios.aroundRatioAll.size(); ++i)
	{
		double value = ratios.aroundRatioAll[i][scale];
			grid[i].proximity_ = value;
			grid[i].multiPurposeValue_ = value;
			allSum += value;

			if (currentDec != grid[i].zenith_){
				currentDec = grid[i].zenith_;
				radialSymmetricValues.push_back(0.0);
				counters.push_back(0);
				//std::cerr << "newly filling in " << radialSymmetricValues.size()-1 << " " << i << std::endl;
			}
			unsigned int index = radialSymmetricValues.size()-1;
			counters[index] += 1;
			radialSymmetricValues[index] += value;
	}

	if (radialSymmetricValues.size() != counters.size()){
		std::cerr << "sizes do not match" << radialSymmetricValues.size() << " " << counters.size() << std::endl;
		throw -3;
	}

	for (unsigned int i = 0; i < radialSymmetricValues.size(); ++i)
	{
		radialSymmetricValues[i] /= static_cast<double>(counters[i]);
	}
}

void smoothRadialSymmetricValues(std::vector<double>& radialSymmetricValues, unsigned int n = 0, std::string outputFolder = "")
{
	if (outputFolder != "") dumpDebugOutputOnZeroOrTwenty(n, radialSymmetricValues, std::string(outputFolder+"radialSymmetricNormalization") );

	// smoothing away the last bins on both ends (don't trust the poles with their small statistic!)
	unsigned int numSmoothAways = 3;
	for (unsigned int i = 0; i < numSmoothAways; ++i)
	{
		radialSymmetricValues[i] = radialSymmetricValues[numSmoothAways];
		radialSymmetricValues[radialSymmetricValues.size()-1-i] = radialSymmetricValues[radialSymmetricValues.size()-1-numSmoothAways];
	}

	if (outputFolder != "") dumpDebugOutputOnZeroOrTwenty(n, radialSymmetricValues, std::string(outputFolder + "radialSymmetricNormalization1Smoothed") );

	sgcu::smoothMedian(radialSymmetricValues);	// this used to be median5 !
	// sgcu::smoothMedian5(radialSymmetricValues);

	if (outputFolder != "")	dumpDebugOutputOnZeroOrTwenty(n, radialSymmetricValues, std::string(outputFolder + "radialSymmetricNormalization2Smoothed") );

	unsigned int numSmoothNoCloseToZeros = 5;	// for ANTARES 2 and IC40 2
//	numSmoothNoCloseToZeros = 3;	// for IceCube, large dataset

	for (unsigned int i = 0; i < numSmoothNoCloseToZeros; ++i)
	{
		smoothNoCloseToZerosKillingEnds(radialSymmetricValues);
	}

	if (outputFolder != "") dumpDebugOutputOnZeroOrTwenty(n, radialSymmetricValues, std::string(outputFolder+"radialSymmetricNormalization3Smoothed") );

	// #CODETODETERMINETHENUMBEROFLOWPASS
	unsigned int numSmoothNoZeros = 60;	// for ANTARES
	numSmoothNoZeros = 15;	// for IceCube, IC40
	numSmoothNoZeros = 2;	// for ANTARES 2 and IC40 2
//	numSmoothNoZeros = 1;	// for IceCube, large dataset
	// std::cout << "debug: smoothRadialSymmetricValues: number of smoothing operations : " << numSmoothNoZeros << std::endl;
	for (unsigned int i = 0; i < numSmoothNoZeros; ++i)
	{
		// sgcu::gaussianSmooth(radialSymmetricValues);	// used to be smoothNoZeros (due to bug without killing ends)
		smoothNoZerosKillingEnds(radialSymmetricValues);
	}

	//sgcu::smooth(radialSymmetricValues);

	if (outputFolder != "") dumpDebugOutputOnZeroOrTwenty(n, radialSymmetricValues, std::string(outputFolder + "radialSymmetricNormalization5Smoothed") );
}

void distributeSymmetricValuesOnSphere(std::vector<double> const & radialSymmetricValues, SG::sphericalGrid<double>& grid)
{
	unsigned int index = 0;
	double currentDec = grid[0].zenith_;
	for (unsigned int i = 0; i < grid.size(); ++i)
	{
		if (grid[i].zenith_ != currentDec)
		{
			currentDec = grid[i].zenith_;
			++index;
		}
		grid[i].proximity_ = radialSymmetricValues[index];
	}
}

void signalFirstInternal::computeOneNormalization(const std::vector<SG::candidate_t>& events, unsigned int fine)
{
	std::string folderBkp = outputFolder_;
	outputFolder_ = "intermediate/step1/";

	timeval startt, endt;
	gettimeofday(&startt, 0);

	// Set up the grid to search and evaluate on
	grid_.setupGrid(fine);
	ratios_.init(stepDegree_);

	// get randomized events
	std::vector<SG::candidate_t> eventsRandomized;
	getRandomEventsFromEvents(events, eventsRandomized);
	myVtkWriter_.write(eventsRandomized, std::string(outputFolder_+"eventsNorm3D") );

	gettimeofday(&endt, 0); std::cout << "Runtime of setup: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);

	// evaluate on every grid point the distance to all events
	ratios_.computeAllRatios(grid_, eventsRandomized);

	gettimeofday(&endt, 0); std::cout << "Runtime of computeAllRatios: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);

	unsigned int numDistanceBins = ratios_.aroundRatioAll[0].size();
	for (unsigned int n = 0; n < numDistanceBins; ++n)
	{
		// Check if the ratios and the grid match in size
		if (ratios_.aroundRatioAll.size() != grid_.size()){
			std::cerr << "Size mismatch ratios: " << ratios_.aroundRatioAll.size() << " vs grid:" << grid_.size() << std::endl;
			return;
		}

		double allSum = 0.0;
		std::vector<double> radialSymmetricValues;
		computeRadialSymmetricValuesOneScale(ratios_, grid_, n, radialSymmetricValues, allSum);

		smoothRadialSymmetricValues(radialSymmetricValues, n, outputFolder_);

		distributeSymmetricValuesOnSphere(radialSymmetricValues, grid_);

		renormalizeGrid(grid_, allSum);

		// store the result for this scale in the full set of ratios
		for (unsigned int i = 0; i < grid_.size(); ++i)
		{
			ratios_.aroundRatioAll[i][n] = grid_[i].proximity_;
		}
	}

	gettimeofday(&endt, 0); std::cout << "Runtime of the computing expectations from ratios: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);

	outputFolder_ = folderBkp;
}

void signalFirstInternal::computeNormalizationsFromDataStep1(const std::vector<SG::candidate_t>& events, unsigned int fine)
{
	std::cout << "Normalization Step1 with " << events.size() << " events" << std::endl;

	std::string outputPath = "intermediate/step1/";
	setOutputFolder(outputPath);
	computeOneNormalization(events, fine);

//	writeNormalizationGrids(folderPrefix+"gridNormalizationNew");
	ratios_.dumpRatiosAsNormalization(outputPath+"normalizationDataStep1", events.size());
}

void signalFirstInternal::init(unsigned int fine)
{
	grid_.setupGrid(fine);
	segmentationPoints_.clear();

	for (unsigned int i = 1; i < 5; ++i)
	{
		segmentationPoints_.push_back( std::pow(0.5, i) );
	}
}

void signalFirstInternal::performEvaluate(std::vector<SG::candidate_t>& events, unsigned int fine) {
	std::string outputBkp = outputFolder_;
	outputFolder_ = "intermediate/step5/";

	std::cout << "Evaluating " << events.size() << " events" << std::endl;

	myVtkWriter_.write(events, outputFolder_+"events3D");
	writeEquatorialProjection(events, outputFolder_+"hammerProjEventsEquatorial.txt");
	writeGalacticProjection(events, outputFolder_+"hammerProjEventsGalactic.txt");

	evaluateCore(events, fine, "");

	printMaxRelevanceSegmentationwise(finalClustersMulti_, grid_, segmentationPoints_);

	std::string folderPrefixStep2 = "intermediate/step2/";
	renormalizeStep2MaximaOnly(finalClustersMulti_, segmentationPoints_, folderPrefixStep2+"normalizationDataStep2", "clusters_segment");

 	std::cout << "Maximal relevance, all segmentations, all clusters, excl. trial factor: " << getMaxRelevance(finalClustersMulti_) << std::endl;
 	printMaxRelevanceSegmentationwise(finalClustersMulti_, grid_, segmentationPoints_);

//	for (unsigned int i = 0; i < finalClustersMulti_.size(); ++i)
//	{
//		drawClusters(finalClustersMulti_[i], grid_);
//		myVtkWriter_.write(grid_.gridpoints_, outputFolder_+"grid3DIntermed" + boost::lexical_cast<std::string>(++globalOutputCounter));
//	}
	// output the stencil of the large cluster
//	drawClusters(finalClustersMulti_[0], grid_);
	//grid_.save("stencilGrid.grid");

//	for (unsigned int i = 0; i < finalClustersMulti_.size(); ++i){
//		std::cerr << "DEBUG drawing debug output" << i << std::endl;
//		drawClusters(finalClustersMulti_[i], grid_);
//		std::cerr << "DEBUG writing debug output" << i << std::endl;
//		myVtkWriter_.write(grid_.gridpoints_, outputFolder_+"grid3DIntermed" + boost::lexical_cast<std::string>(++globalOutputCounter));
//	}

//	std::cerr << "DEBUG renorm3" << std::endl;
	renormalizeClusterRelevances(finalClustersMulti_);

	//std::cerr << "DEBUG finished renorm3" << std::endl;
 	std::cout << "Maximal relevance, all segmentations, all clusters, incl. trial factor, with metrics mask: " << getMaxRelevance(finalClustersMulti_) << std::endl;
 	printMaxRelevanceSegmentationwise(finalClustersMulti_, grid_, segmentationPoints_);
	//std::cerr << "printing clusters in ASCII files" << std::endl;
	printClustersToASCII(finalClustersMulti_, grid_);

	printEventInfoPerCluster(events, grid_, finalClustersMulti_, outputFolder_);

	std::cerr << "Drawing finalized clusters" << std::endl;
	for (unsigned int i = 0; i < finalClustersMulti_.size(); ++i)
	{
//		std::cerr << "Drawing finalized cluster " << i << " / " << finalClustersMulti_.size() << std::endl;
		drawClusters(finalClustersMulti_[i], grid_);
//		myVtkWriter_.write(grid_.gridpoints_, outputFolder_+"grid3DIntermed" + boost::lexical_cast<std::string>(++globalOutputCounter));
////	grid_.writeAsGalacticHammerProj(outputFolder_);
		grid_.writeAsEquatorialHammerProj(outputFolder_);
	}

	outputFolder_ = outputBkp;
}

void setFirstScaleToZero(allScales & scales)
{
	for (unsigned int i = 0; i < scales.aroundRatioAll.size(); ++i)
	{
		scales.aroundRatioAll[i][0] = 0.0;
	}
}

void evaluateRatiosToGrid(SG::allScales const & scales, SG::sphericalGrid<double> & grid, std::string outputFolder)
{
	//TODO: check if grid is actually set up at this stage
	sumScalesToGrid(scales, grid);
	grid.bkpProximity();

	std::cout << "OF: " << outputFolder << std::endl;
//	myVtkWriter_.write(grid.gridpoints_, outputFolder+"grid3DIntermed" + boost::lexical_cast<std::string>(++globalOutputCounter));
	grid.writeAsEquatorialHammerProj(outputFolder);
//	writeGridAsGalacticHammerProj(grid);

	// dump grid without hammerProjection
	save(grid, outputFolder+"/gridComputationCoordsDetailed.grid");
}

void checkConsistentSize(std::vector<std::vector<std::vector<SG::cluster> > > const & allEntries, std::vector<double> const & segmentationPoints)
{
	unsigned int numberOfPseudoExperimentsPerSegmentation = allEntries[0].size();
	for (unsigned int n = 1; n < allEntries.size(); ++n)
	{
		if (allEntries[n].size() != numberOfPseudoExperimentsPerSegmentation)
		{
			std::cerr << "Not all Segmentation have the same number of pseudoexperiments!" <<
					"I expected " << numberOfPseudoExperimentsPerSegmentation << " but for segmentation " << n  <<
					" with threshold " << segmentationPoints[n] << " I got " << allEntries[n].size() <<
					"! Cannot proceed to compute normalization3 from normalization 2 values." << std::endl;
			throw -1;
		}
	}
	if (numberOfPseudoExperimentsPerSegmentation == 0)
	{
		std::cerr << "Failed to load any clusters from the pseudo experiments." << std::endl;
		throw -1;
	}
}

// renormalize all clusters of this one pseudo experiment with all their metrics, overwrite the raw metrics values by the renormalized results
// and return the highest renormalized relevance of any cluster and any metric of this pseudo experiment
double renormalizeOnePseudoExperimentAllClustersAllMetrics(std::vector<std::vector<struct rememberBoostHistogram > > const & histsAllSegmentationsAllRelevances,
		std::vector<double> const & segmentationPoints, unsigned int pseudoexperimentIndex,
		std::vector<std::vector<std::vector<SG::cluster> > > & allClustersFromPseudoExperiments)
{
	std::vector<double> relevances;
	// check all segmentations for one pseudo experiment
	for (unsigned int n = 0; n < segmentationPoints.size(); ++n)
	{
		// std::cerr << "in seg " << n << " / " << segmentationPoints_.size() << std::endl;

		std::vector<SG::cluster> comparisonClusters;
		// Now consider all clusters that have been produced in one pseudo experiment with one segmentation threshold
		comparisonClusters = allClustersFromPseudoExperiments[n][pseudoexperimentIndex];
		if (comparisonClusters.size() == 0)
		{
			std::cerr << "No cluster loaded from segmentation " << n << " at " << segmentationPoints[n] << std::endl;
			continue;
		}

		// now translate the relevance values (size, highest value, etc.) to
		// a p-value (how likely is it to observe such a value or greater) and therefore
		// to a significance (sometimes also called renormalized relevance or pre-trial significance)
		renormalizeCore(comparisonClusters, histsAllSegmentationsAllRelevances[n]);	// renormalize all metrics of all clusters of this pseudo experiment of this segmentation

		// use the renormalized clusters from now on
		allClustersFromPseudoExperiments[n][pseudoexperimentIndex] = comparisonClusters;

		for (unsigned int j = 0; j < comparisonClusters.size(); ++j)
		{	// just to make the relevances vector shorter, inprinciple the maximum is not needed here
			relevances.push_back(*std::max_element(comparisonClusters[j].relevances.begin(), comparisonClusters[j].relevances.end()));
		}
	}
	// store the maximal renormalized relevance value of any cluster and any metric for this pseudo experiment
	return ( *max_element(relevances.begin(), relevances.end()) );
}

//void signalFirstInternal::extendNormalizationFromStep2ToStep3()
void signalFirstInternal::renormalizeAllPseudoexperimentsAllClustersAllMetrics()
{
	std::string inputFolder = "intermediate/step2/normalizationDataStep2";
	std::string searchString1 = "clusters_segment";

	std::cout << "Computing normalization3 distribution from normalization2 pseudo-experiments." << std::endl;
	std::cout << "Loading comparison values ..." << std::endl;

	// Load all comparison events
	std::vector<std::vector<struct rememberBoostHistogram > > histsAllSegmentationsAllRelevances;
	std::vector<std::vector<std::vector<SG::cluster> > > allEntries;
	getComparisonValuesAndHistograms(inputFolder, searchString1, segmentationPoints_, allEntries, histsAllSegmentationsAllRelevances);
	checkConsistentSize(allEntries, segmentationPoints_);

	std::vector<double> maximalValues;
	for (unsigned int i = 0; i < allEntries[0].size(); ++i){	// for all pseudo experiments of this segmentation
		// std::cerr << "Processing pseudo-experiment " << i << " / " << allEntries[0].size() << std::endl;
		maximalValues.push_back( renormalizeOnePseudoExperimentAllClustersAllMetrics(histsAllSegmentationsAllRelevances, segmentationPoints_, i, allEntries) );
	}

	const std::string normalizationFolderOutput = "intermediate/step4/normalizationDataStep3/";
	sgcu::saveVectorToFile(maximalValues, normalizationFolderOutput + "normalizationStep3_Maxima.txt");
}

void signalFirstInternal::evaluateEventsToRatios(const std::vector<SG::candidate_t>& events, unsigned int fine) {
	std::cout << "Evaluate events to ratios" << std::endl;

	timeval startt, endt;
	gettimeofday(&startt, 0);

	// Set up the grid to search and evaluate on
	grid_.setupGrid(fine);
	ratios_.init(stepDegree_);

	gettimeofday(&endt, 0); std::cout << "runtime of init: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);

	// The evaluateDensity sets the ratios, not the grid
	evaluateDensity(events);
	ratios_.assertNoNansAndInfs();

	gettimeofday(&endt, 0); std::cout << "runtime of evaluateDensity: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);

	debugWriteRatios();

	segmentAllRatiosScalewise(ratios_);
	ratios_.assertNoNansAndInfs();

	debugWriteRatios();

	gettimeofday(&endt, 0); std::cout << "runtime of segmentAllRatiosScalewise: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);

	remapAllRatiosNew();
	ratios_.assertNoNansAndInfs();

	gettimeofday(&endt, 0); std::cout << "runtime of remapAllRatiosNew: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);
}

std::vector<cluster> evaluateGridToClusters(SG::sphericalGrid<double> & grid, double keepFraction)
{
	std::cout << "grid non zero before keepTop  " << grid.countNumberOfNonZeroValues() << std::endl;
	grid.keepTop(keepFraction);
	std::cout << "grid non zero after keepTop  " << grid.countNumberOfNonZeroValues() << std::endl;
//	if (doDebugOutput) myVtkWriter_.write(grid.gridpoints_, outputFolder_+"grid3DIntermed" + boost::lexical_cast<std::string>(++globalOutputCounter)); //21

	grid.medianFilter();
//	if (doDebugOutput) myVtkWriter_.write(grid.gridpoints_, outputFolder_+"grid3DIntermed" + boost::lexical_cast<std::string>(++globalOutputCounter)); //21

	return clustersFromGrid(grid);
}

void evaluateGridToSegmentsToClusters(SG::sphericalGrid<double> const & grid, std::vector<double> const & segmentationPoints,
		std::string outputFolder, std::vector<std::vector<cluster> > & finalClustersMulti)
{
	std::cout << "Evaluate grid to segments to clusters" << std::endl;
	grid.assertNoNansAndInfs();

	finalClustersMulti.clear();
	SG::sphericalGrid<double> gridToWork = grid;

	unsigned int thisTimeIndex = time(NULL);
	std::cout << "segmentationPoints_.size() " << segmentationPoints.size() << std::endl;
	for (unsigned int i = 0; i < segmentationPoints.size(); ++i)
	{
		std::cout << i+1 << "/" << segmentationPoints.size() << "  keepFraction=" << segmentationPoints[i] << std::endl;

		gridToWork = grid;
		std::vector<cluster> finalClusters = evaluateGridToClusters(gridToWork, segmentationPoints[i]);

		if (outputFolder != "")
		{
			writeClusters(outputFolder + "/clusters_segment" + boost::lexical_cast<std::string>(segmentationPoints[i])
					+ "_" + boost::lexical_cast<std::string>(thisTimeIndex) + ".txt", finalClusters);
		}

		finalClustersMulti.push_back(finalClusters);
	}
}

void signalFirstInternal::evaluateEventsToSkymap(const std::vector<SG::candidate_t>& events, unsigned int fine){
//	myVtkWriter_.write(events, outputFolder_+"events3D");
//	writeGalacticProjection(events, outputFolder_+"hammerProjEventsGalactic.txt");
	writeEquatorialProjection(events, outputFolder_+"hammerProjEventsEquatorial.txt");
	dumpEvents(events, outputFolder_+"events.txt");

	timeval startt, endt;
	gettimeofday(&startt, 0);

	evaluateEventsToRatios(events, fine);

	gettimeofday(&endt, 0); std::cout << "runtime of evalEventsToRatios: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);

	evaluateRatiosToGrid(ratios_, grid_, outputFolder_);

	writeEquatorialProjection(grid_.gridpoints_, outputFolder_+"hammerProjGridDebugEquatorial.txt");

	gettimeofday(&endt, 0); std::cout << "runtime of evaluateRatiosToGrid: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);
}

void signalFirstInternal::evaluateCore(const std::vector<SG::candidate_t>& events, unsigned int fine, std::string outputFolder)
{
	evaluateEventsToSkymap(events, fine);
	evaluateGridToSegmentsToClusters(grid_, segmentationPoints_, outputFolder, finalClustersMulti_);
}

void signalFirstInternal::evaluateDensity(const std::vector<SG::candidate_t>& events)
{
	setToZero(ratios_);
	unsigned int numBins = static_cast<unsigned int>(maximumDistanceToEvaluateInDegree/stepDegree_)+1;
	if (normalization_.size() == 0)
	{
		std::cout << "Automatically trying to load normalization from step1." << std::endl;
		int ret = normalization_.loadNormalization2(grid_.size(), numBins, events.size(), "intermediate/step1/normalizationDataStep1");
		if (ret != 0)
		{
			std::cerr << "Failed to load normalization for step1 from intermediate/step1/normalizationDataStep1 with error code " << ret << std::endl;
			throw -1;
		}
	}
	if (normalization_.size() == 0)
	{
		std::cerr << "evaluateDensity is only possible if a normalization has been computed and loaded before!" << std::endl;
		throw -2;
	}

	ratios_.computeAllRatios(grid_, events);

	std::cout << countNumberOfNonZeros(ratios_) << " non-zero ratio values" << std::endl;

//	debugWriteRatios();

	ratios_.computePoissonProbabilityForAllRatios(normalization_, events.size());

	invertRatiosAndTakeLog10(ratios_);

	washOutWithinScales(ratios_, grid_);
	washOutBetweenScales(ratios_);
}

void addRandomEventsFromEvents(const std::vector<SG::candidate_t>& events, std::vector<SG::candidate_t>& eventsRandomized)
{
	// use events to generate a lot of pseudo random events here!
	for (unsigned int i = 0; i < events.size(); ++i)
	{
		candidate_t c = events[i];
		c.azimuth_ = sgcu::goodRand(2.0*M_PI);
		eventsRandomized.push_back(c);
	}
}

void getRandomEventsFromEvents(const std::vector<SG::candidate_t>& events, std::vector<SG::candidate_t>& eventsRandomized)
{
	eventsRandomized.clear();
	addRandomEventsFromEvents(events, eventsRandomized);
}

void signalFirstInternal::computeNormalizationsFromDataStep2(const std::vector<SG::candidate_t>& eventsRaw, unsigned int fine, unsigned int numberOfRepetitions){
	std::cout << "Normalization Step2 with " << eventsRaw.size() << " events and " << numberOfRepetitions << " repetitions" << std::endl;

	std::string folderPrefix = "intermediate/step2/";
	setOutputFolder(folderPrefix);
	std::string normalizationFolder = folderPrefix+"normalizationDataStep2";

	for (unsigned int n = 0; n < numberOfRepetitions; ++n){
		std::cout << "Normalizations " << n+1 << "/" << numberOfRepetitions << std::endl;

		std::vector<SG::candidate_t> eventsRandomized;
		getRandomEventsFromEvents(eventsRaw, eventsRandomized);	// only give a new random azimuth

		evaluateCore(eventsRandomized, fine, normalizationFolder);
	}
}

//// This outputs a series of grids , each with the ratios found for one distance bin
//void signalFirstInternal::dumpRatiosForDebug(){
//	std::vector<std::vector<float> > temp;
//	temp.resize(ratios_.aroundRatioAll.size());
//	unsigned int numFirst = 6;
//	unsigned int numSecond = 14;
//	for (unsigned int i = 0; i < temp.size(); ++i){
//		temp[i].resize(numFirst+numSecond+1);
//	}
//	for (unsigned int i = 0; i < temp.size(); ++i){
//		for (unsigned int j = 0; j < numFirst; ++j){
//			temp[i][j] = ratios_.aroundRatioAll[i][j];
//		}
//	}
//	for (unsigned int i = 0; i < temp.size(); ++i){
//		for (unsigned int j = numFirst; j < temp[i].size(); ++j){
//			temp[i][j] = ratios_.aroundRatioAll[i][(j-numFirst)*10];
//		}
//	}
//	for (unsigned int i = 0; i < temp.size(); ++i){
//		unsigned int j = numFirst+numSecond;
//		temp[i][j] = ratios_.aroundRatioAll[i][ratios_.aroundRatioAll[i].size()-1];
//	}
//	myVtkWriter_.write(grid_.gridpoints_, temp, "grid3D");
//}

void washOutWithinScales(allScales & scales, SG::sphericalGrid<double> const & grid)
{
	const unsigned int numberOfGeneralSmearings = 1;
	const unsigned int numSmearCentered = 1;

	SG::sphericalGrid<double> tempGrid = grid;
	for (unsigned int k = 0; k < numberOfGeneralSmearings; ++k)
	{
		// apply spatial filters here for every distance bin
		for (unsigned int i = 0; i < scales.aroundRatioAll[0].size(); ++i)
		{

			// copy ratios to the grid	// TODO check if this can't use a function that can be reused
			for (unsigned int j = 0; j < tempGrid.size(); ++j)
			{
				tempGrid[j].proximity_ = scales.aroundRatioAll[j][i];
			}

			// perform some filters
			for (unsigned int rep = 0; rep < numSmearCentered; ++rep)
				tempGrid.smearCentered();

			// copy values back from the grid to the ratios
			for (unsigned int j = 0; j < tempGrid.size(); ++j)
			{
				scales.aroundRatioAll[j][i] = tempGrid[j].proximity_;
			}
		}
	}
}

void washOutBetweenScales(allScales & scales){	// TODO: exchange tempRatios and scales (saves a copy of the whole thing)
	const unsigned int numDistanceGaussians = 1;

	allScales tempRatios = scales;
	// smooth the results over neighbouring distance bins
	unsigned int numBins = scales.aroundRatioAll[0].size();
	for (unsigned int rep = 0; rep < numDistanceGaussians; ++rep){
		for (unsigned int i = 0; i < scales.aroundRatioAll.size(); ++i){
			// narrow centered smear
			for (unsigned int j = 1; j < numBins-1; ++j){
				tempRatios.aroundRatioAll[i][j] = (scales.aroundRatioAll[i][j-1] + 4.0*scales.aroundRatioAll[i][j] + scales.aroundRatioAll[i][j+1]) / 6.0;
			}
			// approximately correct the first and last
			tempRatios.aroundRatioAll[i][0] = (2.0*scales.aroundRatioAll[i][0] + scales.aroundRatioAll[i][1]) / 3.0;
			tempRatios.aroundRatioAll[i][numBins-1] = (scales.aroundRatioAll[i][numBins-2] + 2.0*scales.aroundRatioAll[i][numBins-1])/3.0;
		}
	}
	scales = tempRatios;
}

std::vector<double> signalFirstInternal::computeRemapGridDistanceFractions()
{
	// std::cout << "Computing the remapping grid numbers how many points are at a certain distance" << std::endl;
	// This is the mapping for every distance how many other gridpoints are at this distance
	std::vector<double> mapping;
	unsigned int numBins = 180;
	mapping.resize(numBins);

	// we compute a mapping for every distance
	for (unsigned int i = 0; i < numBins; ++i)
	{
		// how many gridpoints should be at that distance?
		double tempAngle = static_cast<double>( M_PI*i / (18*grid_.fine()) );
		double densityFactor = sin(tempAngle);
		unsigned int howManyGridpoints = static_cast<unsigned int>(grid_.fine()*36*densityFactor + 0.5);
		if (howManyGridpoints == 0) howManyGridpoints = 1;
		mapping[i] = 1.0/static_cast<double>(howManyGridpoints);
	}
	return mapping;
}

void signalFirstInternal::remapAllRatiosNew(){
//	cout << "why aren't you using remapAllRatiosNewAgain() ???" << std::endl;

	unsigned int numGridpoints = grid_.size();
	if (numGridpoints == 0){
		std::cerr << "Remapping ratios is only possible if there is a grid!" << std::endl;
		return;
	}
	unsigned int numBins = ratios_.aroundRatioAll[0].size();
	std::cout << "Remapping ratios" << std::endl;

	double stepRad = 0.5/180*M_PI;
	allScales tempRatios = ratios_;
	for (unsigned int i = 0; i < numGridpoints; ++i){
		for (unsigned int j = 0; j < numBins; ++j){
			tempRatios.aroundRatioAll[i][j] = 0.0;
		}
	}

	static std::vector<double> distanceFractions = computeRemapGridDistanceFractions();

	// remap proximities on every gridpoint
//	omp_set_num_threads(16);
	#pragma omp parallel for
	for (int i = 0; i < (int)numGridpoints; ++i){
		// check with every other gridpoint if it has exactly the distance the found cluster
		for (unsigned int j = 0; j < numGridpoints; ++j){
			double dist = distanceRad(grid_[i], grid_[j]);

			// what is this distance in gridpoint bins?
			unsigned int binIndex = static_cast<unsigned int>(dist/stepRad+0.5);
			if (binIndex < numBins){
				double proximity = ratios_.aroundRatioAll[i][binIndex];
				// do not consider gridpoints which have NO proximity to redistribute!
				if (proximity <= 0.0){
					continue;
				}
				// else: if they have a proximity (meaning e.g. they are part of a found cluster) then process this gridpoint

				// always set the minimal multipurpose value that helped create the resulting remapped gridpoint
//				double multiPurposeValue = binIndex;

				double proximitySplit = proximity*distanceFractions[binIndex];

				#pragma omp critical
				{
					// then add some of the proximity of this cluster point to the found point in the right distance
					tempRatios.aroundRatioAll[j][binIndex] += proximitySplit;
				}
			}
		}
	}
	// finally store our work
	ratios_ = tempRatios;
}

// TODO speed up the remapping, but nothing really worked so far
//void signalFirstInternal::remapAllRatiosNewAgain(){
//	unsigned int numGridpoints = grid_.size();
//	if (numGridpoints == 0){
//		std::cerr << "Remapping ratios is only possible if there is a grid!" << std::endl;
//		return;
//	}
//	unsigned int numBins = ratios_.aroundRatioAll[0].size();
//	std::cout << "Remapping ratios" << std::endl;
//
//	timeval startt, endt;
//	gettimeofday(&startt, 0);
//
//	double stepRad = ratios_.stepSize_;
//	std::cout << "debug: stepRad = " << stepRad << " , used to be " << 0.5/180.0*M_PI << std::endl;
//	allScales tempRatios = ratios_;
//	for (unsigned int i = 0; i < numGridpoints; ++i){
//		for (unsigned int j = 0; j < numBins; ++j){
//			tempRatios.aroundRatioAll[i][j] = 0.0;
//		}
//	}
//
//	gettimeofday(&endt, 0); std::cout << "runtime of zeroing: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);
//
//	static std::vector<double> distanceFractions = computeRemapGridDistanceFractions();
//
//	gettimeofday(&endt, 0); std::cout << "runtime of computeRemapGridDistanceFractions: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);
//
//	// remap proximities on every gridpoint
////	omp_set_num_threads(16);
//	#pragma omp parallel for
//	for (int i = 0; i < 1000; ++i){
////	for (int i = 0; i < (int)numGridpoints; ++i){
//		// check with every other gridpoint if it has exactly the distance the found cluster
//		for (unsigned int j = 0; j < numGridpoints; ++j){
//			double dist = distanceRad(grid_[i], grid_[j]);
//
//			// what is this distance in gridpoint bins?
//			unsigned int binIndex = static_cast<unsigned int>(dist/stepRad+0.5);
//			if (binIndex < numBins){
//				double proximity = ratios_.aroundRatioAll[i][binIndex];
//				// do not consider gridpoints which have NO proximity to redistribute!
//				if (proximity <= 0.0){
//					continue;
//				}
//				// else: if they have a proximity (meaning e.g. they are part of a found cluster) then process this gridpoint
//
//				// always set the minimal multipurpose value that helped create the resulting remapped gridpoint
////				double multiPurposeValue = binIndex;
//
//				double proximitySplit = proximity*distanceFractions[binIndex];
//
//				#pragma omp critical
//				{
//					// then add some of the proximity of this cluster point to the found point in the right distance
//					tempRatios.aroundRatioAll[j][binIndex] += proximitySplit;
//				}
//			}
//		}
//	}
//
//	gettimeofday(&endt, 0); std::cout << "runtime of remapping core: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);
//
//	// finally store our work
//	ratios_ = tempRatios;
//
//	gettimeofday(&endt, 0); std::cout << "runtime of copying ratios: " << endt.tv_sec - startt.tv_sec << std::endl; gettimeofday(&startt, 0);
//}

// TODO check
void sumScalesToGrid(allScales const & scales, SG::sphericalGrid<double> & grid)
{
	for (unsigned int i = 0; i < scales.aroundRatioAll.size(); ++i)
	{
		grid[i].proximity_ = sgcu::accumulate(scales.aroundRatioAll[i]);
	}
}

void invertRatiosAndTakeLog10(allScales & scales)
{
	for (unsigned int i = 0; i < scales.aroundRatioAll.size(); ++i)
	{
		for (unsigned int j = 0; j < scales.aroundRatioAll[i].size(); ++j)
		{
			scales.aroundRatioAll[i][j] = log10( 1.0/scales.aroundRatioAll[i][j] );
		}
	}
}

void writeRatiosToGrid(allScales const & scales, unsigned int distanceIndex, SG::sphericalGrid<double> & grid)
{
	static unsigned int outputCounter = 0;
	++outputCounter;
	for (unsigned int i = 0; i < grid.size(); ++i)
	{
		grid[i].proximity_ = scales.aroundRatioAll[i][distanceIndex];
	}
	writeEquatorialProjection(grid.gridpoints_, "hammerProjEquatorial_distance" +
			boost::lexical_cast<std::string>(distanceIndex) + "_nr"+boost::lexical_cast<std::string>(outputCounter)+".txt");
//	myVtkWriter_.write(grid_.gridpoints_, outputFolder_ + "grid3DIntermed" + boost::lexical_cast<std::string>(outputCounter));
}

void signalFirstInternal::debugWriteRatios()
{
	//std::cerr << " We don't write debugRatios here anymore!" << std::endl;
	return;

	writeRatiosToGrid(ratios_, 0, grid_);
	writeRatiosToGrid(ratios_, 6, grid_);
	writeRatiosToGrid(ratios_, 20, grid_);
	writeRatiosToGrid(ratios_, 100, grid_);
}

//void signalFirstInternal::discardClustersOutsideOfTemplate(double minFraction){
//	//std::cerr << "DEBUG Starting discard cluster based on stencil" << std::endl;
//	// if there is a stencil: discard every cluster that is not in the stencil (at least to a minimum fraction)
//	if (stencilGrid_.gridpoints_.size() != 0){
//		std::cerr << "Beware! Stencil found! Cutting away everything that is not in this stencil!" << std::endl;
//		//for (unsigned int i = 0; i < finalClustersMulti_.size(); ++i){
//		// we operate on finalClusters_ now, no outer loop neccessary
//			//for (unsigned int j = 0; j < finalClustersMulti_[i].size(); ++j){
//			for (unsigned int j = 0; j < finalClusters_.size(); ++j){
//				// check this one cluster
//				//std::vector<unsigned int> clusterpoints = finalClustersMulti_[i][j].gridpoints;
//				std::vector<unsigned int> clusterpoints = finalClusters_[j].gridpoints;
//				//std::cerr << "Testing cluster no " << j << " of size " << clusterpoints.size() << std::endl;
//
//				unsigned int insideCounter = 0;
//				// count how many pixels are inside
//				for (unsigned int k = 0; k < clusterpoints.size(); ++k){
//					unsigned int index = clusterpoints[k];
//					if (stencilGrid_.gridpoints_[index].proximity_ > 0.0){
//						++insideCounter;
//					}
//				}
//				// compute the overlaps-with-stencil-fraction:
//				double fraction = static_cast<double>(insideCounter)/static_cast<double>(clusterpoints.size());
//				//std::cerr << "Cluster has a fraction of " << fraction << std::endl;
//				// discard all gridpoints for this cluster if it doesn't match well enough
//				if (fraction < minFraction){
//					//std::cerr << "Discarding cluster! ###########################" << std::endl;
//					//finalClustersMulti_[i][j].gridpoints.clear();
//					//finalClusters_[j].gridpoints.clear();
//					//for (unsigned int k = 0; k < finalClustersMulti_[i][j].relevance.size(); ++k){
//					//std::cerr << "relevances before: " ;
//					for (unsigned int k = 0; k < finalClusters_[j].relevances.size(); ++k){
//						//std::cerr << finalClusters_[j].relevance[k] << " " << std::endl;
//						//finalClustersMulti_[i][j].relevance[k] = 0.0;
//						finalClusters_[j].relevances[k] = 0.0;
//					}
//					//std::cerr << std::endl;
//				}
//			}
//		//}
//	}
//	//std::cerr << "DEBUG finished discarding cluster" << std::endl;
//}

//void signalFirstInternal::readTemplate(unsigned int fine, std::string fileName, double significanceThreshold) {
//	stencilGrid_.setupGrid(fine);	// shouldn't be neccessary, but I think it wouldn't work without
//	stencilGrid_.load(fileName);
//	//stencilGrid_.correctEquatorialRightAscensionRange();
//
//	for (unsigned int i = 0; i < stencilGrid_.gridpoints_.size(); ++i){
//		if (stencilGrid_.gridpoints_[i].proximity_ < significanceThreshold){
//			stencilGrid_.gridpoints_[i].proximity_ = 0.0;
//		}
//	}
//	stencilGrid_.writeHammerProjection("hammerProjStencilEquatorial.txt");
//}

// put points closely together on the coordinates grid lines
void writeGridLines(){
	std::cerr << "writing grid Lines" << std::endl;
	const double oneDeg = M_PI/180.0;
	const double delta = 0.5*oneDeg;
	const int decStart = -90;
	const int decEnd = 90;
	const int raStart = -180;
	const int raEnd = 180;

	std::vector<SG::candidate<double> > gridLines;
	SG::candidate<double> tempPoint;
	for (int i = decStart; i <= decEnd; i += 10){	// declination from -90 to 90
		double dec = i*oneDeg;
		tempPoint.zenith_ = dec;
		double ra = raStart*oneDeg;
		while (ra <= raEnd*oneDeg){
			tempPoint.azimuth_ = ra;
			tempPoint.proximity_ = ra;
			gridLines.push_back(tempPoint);
			ra += delta;
		}
	}
	std::cerr << "dec finished" << std::endl;
	for (int i = raStart; i <= raEnd; i += 10){ // right ascension from 0 to 360
		double ra = i*oneDeg;
		tempPoint.azimuth_ = ra;
		tempPoint.proximity_ = ra;
		double dec = decStart*oneDeg;
		while (dec <= decEnd*oneDeg){
			tempPoint.zenith_ = dec;
			gridLines.push_back(tempPoint);
			dec += delta;
		}

	}
	//findRange(gridLines);
	std::cerr << "ra finished" << std::endl;
	writeHammerProjection(gridLines, "gridLines.grid");
	std::cerr << "file written" << std::endl;
}

void writeGridLabelsGalactic(){
	std::cerr << "writing grid Labels" << std::endl;
	const double oneDeg = M_PI/180.0;
//	const double delta = 1.0*oneDeg;
	const int decStart = -90;
	const int decEnd = 90;
	const int raStart = -180;
	const int raEnd = 180;
	const int stepSize = 30;

	std::vector<SG::candidate<double> > gridLines;
	SG::candidate<double> tempPoint;
	for (int i = decStart; i <= decEnd; i += 10){	// declination from -90 to 90
		double dec = i*oneDeg;
		double ra = raStart*oneDeg;
		tempPoint.zenith_ = dec;
		tempPoint.azimuth_ = ra;
		tempPoint.proximity_ = i;
		gridLines.push_back(tempPoint);
	}
	std::cerr << "dec finished" << std::endl;
	for (int i = raStart+stepSize; i <= raEnd-stepSize; i += stepSize){ // right ascension from 0 to 360
		double ra = i*oneDeg;
		double dec = 0.5*(decStart+decEnd)*oneDeg;
		tempPoint.azimuth_ = ra;
		tempPoint.zenith_ = dec;
		tempPoint.proximity_ = i;
		gridLines.push_back(tempPoint);
	}
	//findRange(gridLines);
	std::cerr << "ra finished" << std::endl;
//	writeHammerProjection("gridLabelsGal.grid", gridLines);
	writeHammerProjection(gridLines, "gridLabelsGal.grid");
	std::cerr << "file written" << std::endl;
}

void writeGridLabelsEquatorial(){
	std::cerr << "writing grid Labels" << std::endl;
	const double oneDeg = M_PI/180.0;
//	const double delta = 1.0*oneDeg;
	const int decStart = -90;
	const int decEnd = 90;
	const int raStart = -180;
	const int raEnd = 180;
	const int stepSize = 30;

	std::vector<SG::candidate<double> > gridLines;
	SG::candidate<double> tempPoint;
	for (int i = decStart; i <= decEnd; i += 10){	// declination from -90 to 90
		double dec = i*oneDeg;
		double ra = raStart*oneDeg;
		tempPoint.zenith_ = dec;
		tempPoint.azimuth_ = ra;
		tempPoint.proximity_ = i;
		gridLines.push_back(tempPoint);
	}
	std::cerr << "dec finished" << std::endl;
	for (int i = raStart+stepSize; i <= raEnd-stepSize; i += stepSize){ // right ascension from 0 to 360
		double ra = i*oneDeg;
		double dec = 0.5*(decStart+decEnd)*oneDeg;
		tempPoint.azimuth_ = ra;
		tempPoint.proximity_ = i+180;
		tempPoint.zenith_ = dec;
		gridLines.push_back(tempPoint);
	}
	//findRange(gridLines);
	std::cerr << "ra finished" << std::endl;
//	writeHammerProjection("gridLabelsEquat.grid", gridLines);
	writeHammerProjection(gridLines, "gridLabelsEquat.grid");
	std::cerr << "file written" << std::endl;
}

} // end of namespace SG

max_stars_repo_name: sgeisselsoeder/multiscale

max_issues_repo_name: sgeisselsoeder/multiscale

max_forks_repo_name: sgeisselsoeder/multiscale


text:
#ifndef ANTIAIMBOT_DNN_HPP
#define ANTIAIMBOT_DNN_HPP

#include <dlib/dnn.h>

#include "transform.hpp"

using sample_type = output_vector;

using base_network_type =
                    dlib::fc<2,
                    dlib::relu<dlib::fc<6,
                    dlib::relu<dlib::fc<6,
                    dlib::input<sample_type>
                    >>>>>;

using train_network_type = dlib::loss_multiclass_log<base_network_type>;
using test_network_type = dlib::softmax<base_network_type>;

#endif /* ANTIAIMBOT_DNN_HPP */

ext: hpp

max_stars_repo_name: YashasSamaga/sampml

max_issues_repo_name: YashasSamaga/sampml

max_forks_repo_name: YashasSamaga/sampml


text:
#include "test.h"

#include <stdexcept>

#include "baldr/rapidjson_utils.h"
#include <boost/property_tree/ptree.hpp>

#include "baldr/graphid.h"
#include "baldr/graphreader.h"
#include "meili/traffic_segment_matcher.h"

#if !defined(VALHALLA_SOURCE_DIR)
#define VALHALLA_SOURCE_DIR
#endif

using namespace valhalla;

// NOTE: this test relies on pregenerated data. If you change the data format
// you need to update the tiles. The data itself is checked in:
// test/traffic_matcher_tiles/map.osm.gz
// simply gunzip it, create the tiles with it using osmconvert and valhalla_build_tiles
// then individually gzip the tiles using some combination of find xargs and gzip

namespace {

// here we hijack a couple of methods and save off some state while we're at it
// this way the standard calling pattern used from the outside is the same as in the test
// but we now have the internal state so we can see what is going on at more detail
class testable_matcher : public meili::TrafficSegmentMatcher {
public:
  using meili::TrafficSegmentMatcher::TrafficSegmentMatcher;

  std::list<std::vector<meili::interpolation_t>>
  interpolate_matches(const std::vector<meili::MatchResult>& r,
                      std::vector<meili::EdgeSegment>& e,
                      const std::shared_ptr<meili::MapMatcher>& m) const override {
    matches = r;
    matcher = m;
    interpolations = meili::TrafficSegmentMatcher::interpolate_matches(r, e, m);
    return interpolations;
  }

  std::vector<meili::traffic_segment_t>
  form_segments(const std::list<std::vector<meili::interpolation_t>>& i,
                baldr::GraphReader& r) const override {
    segments = meili::TrafficSegmentMatcher::form_segments(i, r);
    return segments;
  }

  mutable std::vector<valhalla::meili::MatchResult> matches;
  mutable std::shared_ptr<meili::MapMatcher> matcher;
  mutable std::list<std::vector<meili::interpolation_t>> interpolations;
  mutable std::vector<meili::traffic_segment_t> segments;
};

// TODO: build the test tiles in the test, need to move traffic association into library to do that
// currently all the logic is in the application

using ots_t = meili::traffic_segment_t; // id,start time,start idx,end time,end idx,length
using ots_matches_t = std::vector<ots_t>;
using sid_t = baldr::GraphId;
std::vector<std::pair<std::string, ots_matches_t>> test_cases{
    // partial, partial with duplicate point
    std::make_pair(
        R"({"trace":[{"lon":-76.376045,"lat":40.539207,"time":0},{"lon":-76.357056,"lat":40.541309,"time":99},{"lon":-76.357056,"lat":40.541309,"time":100}],"trace_options":{"breakage_distance":10000}})",
        ots_matches_t{ots_t{sid_t(0), -1, 0, 50.f, 0, -1}, ots_t{sid_t(0), 50.f, 0, -1, 2, -1}}),
    // partial, partial
    std::make_pair(
        R"({"trace":[{"lon":-76.376045,"lat":40.539207,"time":0},{"lon":-76.357056,"lat":40.541309,"time":100}],"trace_options":{"breakage_distance":10000}})",
        ots_matches_t{ots_t{sid_t(0), -1, 0, 50.f, 0, -1}, ots_t{sid_t(0), 50.f, 0, -1, 1, -1}}),
    // partial, full, partial
    std::make_pair(
        R"({"trace":[{"lon":-76.376045,"lat":40.539207,"time":0},{"lon":-76.351089,"lat":40.541504,"time":300}],"trace_options":{"breakage_distance":10000}})",
        ots_matches_t{ots_t{sid_t(0), -1, 0, 110.f, 0, -1}, ots_t{sid_t(0), 110.f, 0, 250.f, 0, 1000},
                      ots_t{sid_t(0), 250.f, 0, -1, 1, -1}}),
    // partial, full, full, full
    std::make_pair(
        R"({"trace":[{"lon":-76.38126,"lat":40.55602,"time":0},{"lon":-76.35784,"lat":40.56786,"time":600}],"trace_options":{"breakage_distance":10000}})",
        ots_matches_t{ots_t{sid_t(0), -1, 0, 60.f, 0, -1}, ots_t{sid_t(0), 60.f, 0, 110.f, 0, 200},
                      ots_t{sid_t(0), 110.f, 0, 350.f, 0, 1000},
                      ots_t{sid_t(0), 350.f, 0, 600.f, 1, 1000}}),
    // full, full, partial
    std::make_pair(
        R"({"trace":[{"lon":-76.35784,"lat":40.56786,"time":0},{"lon":-76.38126,"lat":40.55602,"time":600}],"trace_options":{"breakage_distance":10000}})",
        ots_matches_t{ots_t{sid_t(0), 0.f, 0, 250.f, 0, 1000},
                      ots_t{sid_t(0), 250.f, 0, 490.f, 0, 1000},
                      ots_t{sid_t(0), 490.f, 0, -1, 0, -1}}),

    // TODO: add test where its all full segments
    // TODO: add test where you are on at the start of a segment, you get off on a small road in
    // between,
    // but come back on again before the segment ends, this segment should be seen twice in output
    // as partials
    // TODO: add test where you are consecutively in the same spot at different times, ie you aren't
    // moving
    // TODO: add test where there is discontinuity in matches so it has to do two sets of matches
    // TODO: add test where intermediate trace points dont get matches, this causes their times to
    // not be used
    // for interpolation but we can still get valid segments on the edges for the entire trace
    // TODO: add a test where you enter a segment, leave it and come back onto it where it starts,
    // via loop,
    // then finish it and you should see partial, then full and the full should not count the length
    // of the partial in it
};

TEST(TrafficMatcher, test_matcher) {
  // fake config
  std::stringstream conf_json;
  conf_json << R"({
      "mjolnir":{"tile_dir":"test/traffic_matcher_tiles"},
      "meili":{"customizable": ["breakage_distance"],
               "mode":"auto","grid":{"cache_size":100240,"size":500},
               "default":{"beta":3,"breakage_distance":2000,"geometry":false,"gps_accuracy":5.0,"interpolation_distance":10,
               "max_route_distance_factor":5,"max_route_time_factor":5,"max_search_radius":100,"route":true,
               "search_radius":50,"sigma_z":4.07,"turn_penalty_factor":200}}
    })";
  boost::property_tree::ptree conf;
  rapidjson::read_json(conf_json, conf);
  conf.get_child("mjolnir").put("tile_dir", VALHALLA_SOURCE_DIR "test/traffic_matcher_tiles");

  // find me a find, catch me a catch
  testable_matcher matcher(conf);

  // some edges should have no matches and most will have no segments
  for (const auto& test_case : test_cases) {
    auto json = matcher.match(test_case.first);
    std::stringstream json_ss;
    json_ss << json;
    boost::property_tree::ptree answer;
    rapidjson::read_json(json_ss, answer);

    const auto& a_segs = test_case.second;
    auto& b_segs = matcher.segments;

    // TODO: for portions of routes with no ots coverage add matching fake segments
    // for now we just remove the fake ones from the result
    auto seg_itr =
        std::remove_if(b_segs.begin(), b_segs.end(), [](const meili::traffic_segment_t& ots) {
          return !ots.segment_id.Is_Valid();
        });
    b_segs.erase(seg_itr, b_segs.end());

    ASSERT_EQ(a_segs.size(), b_segs.size());
    for (size_t i = 0; i < a_segs.size(); ++i) {
      const auto& a = a_segs[i];
      const auto& b = b_segs[i];
      EXPECT_EQ(a.begin_shape_index, b.begin_shape_index);
      EXPECT_EQ(a.end_shape_index, b.end_shape_index);
      EXPECT_NEAR(a.start_time, b.start_time, 10) << "start time is out of tolerance";
      EXPECT_NEAR(a.end_time, b.end_time, 10) << "end time is out of tolerance";
      EXPECT_NEAR(a.length, b.length, 50) << "length is out of tolerance";
    }
  }
}

} // namespace

int main(int argc, char* argv[]) {
  testing::InitGoogleTest(&argc, argv);
  return RUN_ALL_TESTS();
}

ext: cc

max_stars_repo_name: stadiamaps/valhalla

max_issues_repo_name: stadiamaps/valhalla

max_forks_repo_name: eric-erki/valhalla


text:
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements.  See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership.  The ASF licenses this file
// to you under the Apache License, Version 2.0 (the
// "License"); you may not use this file except in compliance
// with the License.  You may obtain a copy of the License at
//
//   http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing,
// software distributed under the License is distributed on an
// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied.  See the License for the
// specific language governing permissions and limitations
// under the License.
//
// The following only applies to changes made to this file as part of YugaByte development.
//
// Portions Copyright (c) YugaByte, Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
// in compliance with the License.  You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the License
// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
// or implied.  See the License for the specific language governing permissions and limitations
// under the License.
//
// Copyright(c) 2020 Futurewei Cloud
//
// Permission is hereby granted,
//        free of charge, to any person obtaining a copy of this software and associated documentation files(the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions :
//
// The above copyright notice and this permission notice shall be included in all copies
// or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS",
// WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
//        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
//        DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//

#include <boost/algorithm/string.hpp>

#include "common/type/decimal.h"
#include "common/k2pg-internal.h"
#include "pggate/pg_op.h"
#include "pggate/pg_env.h"
#include "pggate/pg_gate_typedefs.h"
#include "pggate/catalog/sql_catalog_defaults.h" // for the table/index name constants
#include <string>

namespace k2pg {
namespace gate {
using k2pg::Decimal;

template <typename T>  // static checker to see if a type is numeric (i.e. we can get it out from field by value)
constexpr bool isNumericType() { return std::is_arithmetic<T>::value || std::is_enum<T>::value; }

// template resolution for types that are not supported
template <typename T>
std::enable_if_t<!isNumericType<T>(), Status>
TranslateUserCol(int index, const K2PgTypeEntity* type_entity, const PgTypeAttrs* type_attrs, std::optional<T> field, PgTuple* pg_tuple) {
    return STATUS(InternalError, "unsupported type for user column");
}

// translate numeric types (integers, bool, floats)
template <typename T>
std::enable_if_t<isNumericType<T>(), Status>
TranslateUserCol(int index, const K2PgTypeEntity* type_entity, const PgTypeAttrs* type_attrs, std::optional<T> field, PgTuple* pg_tuple) {
    switch (type_entity->k2pg_type) {
        case K2SQL_DATA_TYPE_INT8: {
            int8_t val = (int8_t)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_INT16: {
            int16_t val = (int16_t)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_INT32: {
            int32_t val = (int32_t)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_INT64: {
            int64_t val = (int64_t)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_UINT32: {
            uint32_t val = (uint32_t)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_UINT64: {
            uint64_t val = (uint64_t)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_BOOL: {
            bool val = (bool)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_FLOAT: {
            float val = (float)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_DOUBLE: {
            double val = (double)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_TIMESTAMP: {
            int64_t val = (int64_t)field.value();
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(&val, sizeof(val), type_attrs));
            break;
        }
        default:
            K2LOG_E(log::pg, "Internal error: unsupported type {}", type_entity->k2pg_type);
            return STATUS(InternalError, "unsupported type for user column");
    }
    return Status::OK();
}

// translate k2::String -based types
template <>
Status
TranslateUserCol<k2::String>(int index, const K2PgTypeEntity* type_entity, const PgTypeAttrs* type_attrs, std::optional<k2::String> field, PgTuple* pg_tuple) {
    switch (type_entity->k2pg_type) {
        case K2SQL_DATA_TYPE_BINARY: {
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(field.value().c_str(), field.value().size(), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_STRING: {
            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(field.value().c_str(), field.value().size(), type_attrs));
            break;
        }
        case K2SQL_DATA_TYPE_DECIMAL: {
            // TODO use SKV/c++ -native decimal64 type
            std::string serialized_decimal(field.value().c_str(), field.value().size());
            Decimal k2pg_decimal;
            if (!k2pg_decimal.DecodeFromComparable(serialized_decimal).ok()) {
                K2LOG_E(log::pg, "Failed to deserialize DECIMAL from {}", serialized_decimal);
                return STATUS(InternalError, "failed to deserialize DECIMAL");
            }
            auto plaintext = k2pg_decimal.ToString();

            pg_tuple->WriteDatum(index, type_entity->k2pg_to_datum(plaintext.c_str(), field.value().size(), type_attrs));
            break;
        }
        default:
            K2LOG_E(log::pg, "Internal error: unsupported type {}", type_entity->k2pg_type);
            return STATUS(InternalError, "unsupported type for user column");
    }
    return Status::OK();
}

template<typename T>
Status TranslateSysCol(int attr_num, std::optional<T> field, PgTuple* pg_tuple) {
    return STATUS(InternalError, "unsupported type for system column");
}

template<>
Status TranslateSysCol<int64_t>(int attr_num, std::optional<int64_t> field, PgTuple* pg_tuple) {
    switch (attr_num) {
        case static_cast<int>(PgSystemAttrNum::kSelfItemPointer):
            pg_tuple->syscols()->ctid = field.value();
            break;
        case static_cast<int>(PgSystemAttrNum::kObjectId):
            pg_tuple->syscols()->oid = field.value();
            break;
        case static_cast<int>(PgSystemAttrNum::kMinTransactionId):
            pg_tuple->syscols()->xmin = field.value();
            break;
        case static_cast<int>(PgSystemAttrNum::kMinCommandId):
            pg_tuple->syscols()->cmin = field.value();
            break;
        case static_cast<int>(PgSystemAttrNum::kMaxTransactionId):
            pg_tuple->syscols()->xmax = field.value();
            break;
        case static_cast<int>(PgSystemAttrNum::kMaxCommandId):
            pg_tuple->syscols()->cmax = field.value();
            break;
        case static_cast<int>(PgSystemAttrNum::kTableOid):
            pg_tuple->syscols()->tableoid = field.value();
            break;
        default:
            return STATUS(InternalError, "system column is not int64_t compatible");
    }
    return Status::OK();
}

template <>
Status TranslateSysCol<k2::String>(int attr_num, std::optional<k2::String> field, PgTuple* pg_tuple) {
    switch (attr_num) {
        case static_cast<int>(PgSystemAttrNum::kPgTupleId): {
            k2::String& val = field.value();
            pg_tuple->Write(&pg_tuple->syscols()->k2pgctid, (const uint8_t*)val.c_str(), val.size());
            break;
        }
        case static_cast<int>(PgSystemAttrNum::kPgIdxBaseTupleId): {
            k2::String& val = field.value();
            pg_tuple->Write(&pg_tuple->syscols()->k2pgbasectid, (const uint8_t*)val.c_str(), val.size());
            break;
        }
        default:
            return STATUS(InternalError, "system column is not string compatible");
    }
    return Status::OK();
}

template<typename T>
void FieldParser(std::optional<T> field, const k2::String& fieldName, const std::unordered_map<std::string, PgExpr*>& targets_by_name, PgTuple* pg_tuple, Status& result, int32_t* num) {
    auto iter = targets_by_name.find(fieldName.c_str());
    if (iter == targets_by_name.end()) {
        if (k2pg::sql::catalog::CatalogConsts::TABLE_ID_COLUMN_NAME == fieldName.c_str() ||
            k2pg::sql::catalog::CatalogConsts::INDEX_ID_COLUMN_NAME == fieldName.c_str()) {
        }
        else {
            K2LOG_D(log::pg, "Encountered field {}, without target reference", fieldName);
        }
        result = Status::OK();
        return;
    }
    if (!iter->second->is_colref()) {
        result= STATUS(InternalError, "Unexpected expression, only column refs supported in SKV");
        return;
    }
    const PgColumnRef* target = (PgColumnRef*)iter->second;
    int attr_num = target->attr_num();

    if (attr_num < 0) {
        if (!field) {
            result= STATUS(InternalError, "Null system column encountered");
            return;
        }
        result = TranslateSysCol(attr_num, std::move(field), pg_tuple);
    }
    else {
        if (!field) {
            pg_tuple->WriteNull(attr_num-1);
            result = Status::OK();
        }
        else {
            result = TranslateUserCol(attr_num-1, target->type_entity(), target->type_attrs(), std::move(field), pg_tuple);
        }
    }
    (*num)++;
    K2LOG_V(log::pg, "Parsed field {}, num: {} in targets_by_name map with size: {}",
            fieldName, (*num), targets_by_name.size());
}

PgOpResult::PgOpResult(std::vector<k2::dto::SKVRecord>&& data) : data_(std::move(data)) {
    ProcessSystemColumns();
}

PgOpResult::PgOpResult(std::vector<k2::dto::SKVRecord>&& data, std::list<int64_t>&& row_orders):
    data_(std::move(data)), row_orders_(move(row_orders)) {
    ProcessSystemColumns();
}

PgOpResult::~PgOpResult() {
}

int64_t PgOpResult::NextRowOrder() {
    return row_orders_.size() > 0 ? row_orders_.front() : -1;
}

// Get the postgres tuple from this batch.
Status PgOpResult::WritePgTuple(const std::vector<PgExpr *> &targets, const std::unordered_map<std::string, PgExpr*>& targets_by_name, PgTuple *pg_tuple, int64_t *row_order) {
    Status result;
    K2ASSERT(log::pg, targets_by_name.size() > 0, "targets should not be empty");
    K2ASSERT(log::pg, syscol_processed_, "System columns have not been processed yet");
    int32_t num = 0;
    FOR_EACH_RECORD_FIELD(data_[nextToConsume_], FieldParser, targets_by_name, pg_tuple, result, &num);
     if (targets_by_name.find("k2pgctid") != targets_by_name.end()) {
        // k2pgctid is a virtual column and won't be in the SKV record
        num++;
    }
    K2ASSERT(log::pg, num == targets_by_name.size(), "All target columns should be processed: {} != {}", num, targets_by_name.size());

    if (pg_tuple->syscols()) {
        auto& k2pgctid_str = k2pgctid_strings_[nextToConsume_];
        pg_tuple->syscols()->k2pgctid = (uint8_t*)k2pg::K2PgCStringToTextWithLen(k2pgctid_str.data(), k2pgctid_str.size());
    }
    K2LOG_D(log::pg, "wrote tuple k2pgctid={}", k2pgctid_strings_[nextToConsume_]);
    if (row_orders_.size()) {
        *row_order = row_orders_.front();
        row_orders_.pop_front();
    } else {
        *row_order = -1;
    }
    ++nextToConsume_;

    return result;
}

// For secondary index read result, where the caller need to get all base row's k2pgctid
void PgOpResult::GetBaseRowIdBatch(std::vector<std::string>& baseRowIds) {
    if (!is_eof())
    {
        for (k2::dto::SKVRecord& record : data_) {
            std::optional<k2::String> basek2pgctid = record.deserializeField<k2::String>("k2pgidxbasectid");

            if (!basek2pgctid.has_value()) {
                CHECK(basek2pgctid.has_value()) << "k2pgidxbasectid for index row was null";
            }
            baseRowIds.emplace_back(*basek2pgctid);
            ++nextToConsume_;
        }
    }
    else
    {
        K2ASSERT(log::pg, false, "BugBug try to get baseRowIdBatch from empty index result!");
    }
}

// Get system columns' values from this batch.
// Currently, we only have k2pgctids, but there could be more.
Status PgOpResult::ProcessSystemColumns() {
    if (syscol_processed_) {
        return Status::OK();
    }
    syscol_processed_ = true;
    for (auto& rec: data_) {
        k2pgctid_strings_.push_back(K2Adapter::GetRowIdFromReadRecord(rec));
        k2pgctids_.emplace_back(k2pgctid_strings_.back().c_str(), k2pgctid_strings_.back().size());
    }
    return Status::OK();
}

//--------------------------------------------------------------------------------------------------

PgOp::PgOp(const std::shared_ptr<PgSession>& pg_session,
                const std::shared_ptr<PgTableDesc>& table_desc,
                const PgObjectId& relation_id)
    : pg_session_(pg_session), table_desc_(table_desc), relation_id_(relation_id) {
    exec_params_.limit_count = 0;
    exec_params_.limit_offset = 0;
    exec_params_.limit_use_default = false;
}

PgOp::~PgOp() {
    // Wait for result in case request was sent.
    // Operation can be part of transaction it is necessary to complete it before transaction commit.
    if (requestAsyncRunResult_.valid()) {
        K2LOG_D(log::pg, "Waiting for in progress request");
        __attribute__((unused)) auto status = requestAsyncRunResult_.get();
    }
}

void PgOp::ExecuteInit(const PgExecParameters *exec_params) {
    end_of_data_ = false;
    if (exec_params) {
        exec_params_ = *exec_params;
    }
}

Result<bool> PgOp::Execute() {
    // SKV is stateless and we have to call query execution every time, i.e., Exec & Fetch, Exec & Fetch
    exec_status_ = SendRequest();
    RETURN_NOT_OK(exec_status_);
    return requestAsyncRunResult_.valid();
}

Status PgOp::GetResult(std::list<PgOpResult> *rowsets) {
    // If the execution has error, return without reading any rows.
    RETURN_NOT_OK(exec_status_);

    if (!end_of_data_) {
        // Send request now in case prefetching was suppressed.
        if (suppress_next_result_prefetching_ && !requestAsyncRunResult_.valid()) {
            K2LOG_D(log::pg, "suppress_next_result_prefetching_: {} send request...", suppress_next_result_prefetching_);
            exec_status_ = SendRequest();
            RETURN_NOT_OK(exec_status_);
        }

        DCHECK(requestAsyncRunResult_.valid());
        auto rows = VERIFY_RESULT(ProcessResponse(requestAsyncRunResult_.get()));
        // In case ProcessResponse doesn't fail with an error
        // it should return non empty rows and/or set end_of_data_.
        DCHECK(!rows.empty() || end_of_data_);
        K2LOG_D(log::pg, "GetResult rows: {}, end_of_data_: {}", rows.size(), end_of_data_);
        rowsets->splice(rowsets->end(), rows);
        // Prefetch next portion of data if needed.
        if (!(end_of_data_ || suppress_next_result_prefetching_)) {
            K2LOG_D(log::pg,
                "Not done, end_of_data_: {}, suppress_next_result_prefetching_:{}, continue sending request...",
                end_of_data_, suppress_next_result_prefetching_);
            exec_status_ = SendRequest();
            RETURN_NOT_OK(exec_status_);
        }
    } else {
        K2LOG_D(log::pg, "Done, GetResult end_of_data_: {}", end_of_data_);
    }

    return Status::OK();
}

Result<int32_t> PgOp::GetRowsAffectedCount() const {
    RETURN_NOT_OK(exec_status_);
    DCHECK(end_of_data_);
    return rows_affected_count_;
}

Status PgOp::ClonePgsqlOps(int op_count) {
    SCHECK(op_count > 0, InternalError, "Table must have at least one partition");
    if (pgsql_ops_.size() < op_count) {
        pgsql_ops_.resize(op_count);
        for (int idx = 0; idx < op_count; idx++) {
        pgsql_ops_[idx] = CloneFromTemplate();

        // Initialize as inactive. Turn it on when setup argument for a specific partition.
        pgsql_ops_[idx]->set_active(false);
        }

        // Set parallism_level_ to maximum possible of operators to be executed at one time.
        parallelism_level_ = pgsql_ops_.size();
    }

    return Status::OK();
}

void PgOp::MoveInactiveOpsOutside() {
    // Move inactive op to the end.
    const int total_op_count = pgsql_ops_.size();
    bool has_sorting_order = !batch_row_orders_.empty();
    int left_iter = 0;
    int right_iter = total_op_count - 1;
    while (true) {
        // Advance left iterator.
        while (left_iter < total_op_count && pgsql_ops_[left_iter]->is_active()) left_iter++;

        // Advance right iterator.
        while (right_iter >= 0 && !pgsql_ops_[right_iter]->is_active()) right_iter--;

        // Move inactive operator to the end by swapping the pointers.
        if (left_iter < right_iter) {
        std::swap(pgsql_ops_[left_iter], pgsql_ops_[right_iter]);
        if (has_sorting_order) {
            std::swap(batch_row_orders_[left_iter], batch_row_orders_[right_iter]);
        }
        } else {
            break;
        }
    }

    // Set active op count.
    active_op_count_ = left_iter;
}

Status PgOp::SendRequest() {
    DCHECK(exec_status_.ok());
    DCHECK(!requestAsyncRunResult_.valid());
    DCHECK(!end_of_data_);
    exec_status_ = SendRequestImpl();
    return exec_status_;
}

Status PgOp::SendRequestImpl() {
    // Populate collected information into requests before sending to SKV.
    RETURN_NOT_OK(CreateRequests());

    // Send at most "parallelism_level_" number of requests at one time.
    int32_t send_count = std::min(parallelism_level_, active_op_count_);
    requestAsyncRunResult_ = VERIFY_RESULT(pg_session_->RunAsync(pgsql_ops_.data(), send_count, relation_id_, &read_time_));

    return Status::OK();
}

Result<std::list<PgOpResult>> PgOp::ProcessResponse(const Status& status) {
    // Check operation status.
    DCHECK(exec_status_.ok());
    exec_status_ = status;
    if (exec_status_.ok()) {
        auto result = ProcessResponseImpl();
        if (result.ok()) {
        return result;
        }
        exec_status_ = result.status();
    }
    return exec_status_;
}

Result<std::list<PgOpResult>> PgOp::ProcessResponseResult() {
    K2LOG_D(log::pg, "Received response for request");

    // Check for errors reported by storage server.
    for (int op_index = 0; op_index < active_op_count_; op_index++) {
        RETURN_NOT_OK(pg_session_->HandleResponse(*pgsql_ops_[op_index], PgObjectId()));
    }

    // Process data coming from storage server.
    std::list<PgOpResult> result;
    bool no_sorting_order = batch_row_orders_.size() == 0;

    rows_affected_count_ = 0;
    for (int op_index = 0; op_index < active_op_count_; op_index++) {
        PgOpTemplate *pgsql_op = pgsql_ops_[op_index].get();
        // Get total number of rows that are operated on.
        rows_affected_count_ += pgsql_op->response().rows_affected_count;

        // Get contents.
        if (!pgsql_op->rows_data().empty()) {
        if (no_sorting_order) {
            result.emplace_back(pgsql_op->rows_data());
        } else {
            result.emplace_back(pgsql_op->rows_data(), std::move(batch_row_orders_[op_index]));
        }
        }
    }

    return result;
}

void PgOp::AddExpr(std::unique_ptr<PgExpr> expr) {
  exprs_.push_back(std::move(expr));
}

//-------------------------------------------------------------------------------------------------

PgReadOp::PgReadOp(const std::shared_ptr<PgSession>& pg_session,
                        const std::shared_ptr<PgTableDesc>& table_desc,
                        std::unique_ptr<PgReadOpTemplate> read_op)
    : PgOp(pg_session, table_desc), template_op_(std::move(read_op)) {
}

void PgReadOp::ExecuteInit(const PgExecParameters *exec_params) {
    PgOp::ExecuteInit(exec_params);

    template_op_->set_return_paging_state(true);
    SetRequestTotalLimit();
    SetRowMark();
    SetReadTime();
}

void PgReadOp::SetReadTime() {
    read_time_ = exec_params_.read_time;
}

Result<std::list<PgOpResult>> PgReadOp::ProcessResponseImpl() {
    // Process result from storage server and check result status.
    auto result = VERIFY_RESULT(ProcessResponseResult());

    // Process paging state and check status.
    RETURN_NOT_OK(ProcessResponsePagingState());
    K2LOG_D(log::pg, "ProcessResponseImpl for ReadOp with result size: {}", result.size());
    return result;
}

Status PgReadOp::CreateRequests() {
    if (request_population_completed_) {
        return Status::OK();
    }

    // No optimization.
    // TODO: create separate requests for different partitions once SKV partition information is available
    pgsql_ops_.push_back(template_op_);
    template_op_->set_active(true);
    active_op_count_ = 1;
    request_population_completed_ = true;
    return Status::OK();
}

Status PgReadOp::InitializeRowIdOperators() {
    // we only support one partition for now
    // keep this logic so that we could support multiple partitions in the future
    int op_count = table_desc_->GetPartitionCount();

    if (batch_row_orders_.size() == 0) {
        // First batch:
        // - Create operators.
        // - Allocate row orders for each tablet server.
        // - Protobuf fields in requests are not yet set so not needed to be cleared.
        RETURN_NOT_OK(ClonePgsqlOps(op_count));
        batch_row_orders_.resize(op_count);
    } else {
        // Second and later batches: Reuse all state variables.
        // - Clear row orders for this batch to be set later.
        // - Clear protobuf fields k2pgctids and others before reusing them in this batch.
        RETURN_NOT_OK(ResetInactivePgsqlOps());
    }
    return Status::OK();
}

Status PgReadOp::PopulateDmlByRowIdOps(const std::vector<std::string>& k2pgctids) {
    // This function is called only when k2pgctids were returned from INDEX, for example,
    //    SELECT xxx FROM <table> WHERE k2pgctid IN (SELECT k2pgctid FROM INDEX);

    RETURN_NOT_OK(InitializeRowIdOperators());
    // Begin a batch of k2pgctids.
    end_of_data_ = false;

    // we only have one partition so far
    PgReadOpTemplate *read_op = GetReadOp(0);
    read_op->set_active(true);
    std::shared_ptr<SqlOpReadRequest> request = read_op->request();
    const K2PgTypeEntity *string_type = K2PgFindTypeEntity(STRING_TYPE_OID);
    // populate k2pgctid values.
    request->k2pgctid_column_values.clear();
    for (const std::string& k2pgctid : k2pgctids) {
        std::unique_ptr<PgConstant> pg_const = std::make_unique<PgConstant>(string_type, SqlValue(k2pgctid));
        // use one batch for now, could split into multiple batches later for optimization
        request->k2pgctid_column_values.push_back(std::make_shared<BindVariable>(static_cast<int32_t>(PgSystemAttrNum::kPgTupleId), pg_const.get()));
        // add to the expr list so that it could be released once PgOP is out of scope
        AddExpr(std::move(pg_const));
    }
    K2LOG_D(log::pg, "Populated {} k2pgctids in op read request for table {}",
            request->k2pgctid_column_values.size(), request->table_id);

    // Done creating request
    MoveInactiveOpsOutside();
    request_population_completed_ = true;

    return Status::OK();
}

Status PgReadOp::ProcessResponsePagingState() {
    K2LOG_D(log::pg, "Processing response paging state for PgReadOp for table {}", relation_id_.GetTableUuid());
    // For each read_op, set up its request for the next batch of data or make it in-active.
    bool has_more_data = false;
    int32_t send_count = std::min(parallelism_level_, active_op_count_);

    for (int op_index = 0; op_index < send_count; op_index++) {
        PgReadOpTemplate *read_op = GetReadOp(op_index);
        auto& res = read_op->response();
        // Check for completion.
        bool has_more_arg = false;
        read_op->request()->paging_state = res.paging_state;
        if (res.paging_state != nullptr) {
            has_more_arg = true;
            K2LOG_D(log::pg, "Updated paging state for request for table {}", read_op->request()->table_id);
        } else {
            K2LOG_D(log::pg, "Response paging state is null for table {}", read_op->request()->table_id);
        }

        if (has_more_arg) {
            has_more_data = true;
            K2LOG_D(log::pg, "has_more_arg -> has_more_data = true");
        } else {
            read_op->set_active(false);
        }
    }

    if (has_more_data || send_count < active_op_count_) {
        // Move inactive ops to the end of pgsql_ops_ to make room for new set of arguments.
        MoveInactiveOpsOutside();
        end_of_data_ = false;
        K2LOG_D(log::pg, "end_of_data = false with has_more_data: {}, send_count: {}, active_op_count: {}",
                has_more_data, send_count, active_op_count_);
    } else {
        // There should be no active op left in queue.
        active_op_count_ = 0;
        end_of_data_ = request_population_completed_;
        K2LOG_D(log::pg, "active_op_count: 0, end_of_data: {}", end_of_data_);
    }

    return Status::OK();
}

void PgReadOp::SetRequestPrefetchLimit() {
    // Predict the maximum prefetch-limit
    std::shared_ptr<SqlOpReadRequest> req = template_op_->request();
    int predicted_limit = default_psql_prefetch_limit;
    if (!req->is_forward_scan) {
        // Backward scan is slower than forward scan, so predicted limit is a smaller number.
        predicted_limit = predicted_limit * default_psql_backward_prefetch_scale_factor;
    }

    // System setting has to be at least 1 while user setting (LIMIT clause) can be anything that
    // is allowed by SQL semantics.
    if (predicted_limit < 1) {
        predicted_limit = 1;
    }

    // Use statement LIMIT(count + offset) if it is smaller than the predicted limit.
    int64_t limit_count = exec_params_.limit_count + exec_params_.limit_offset;
    suppress_next_result_prefetching_ = true;
    K2LOG_D(log::pg, "suppress_next_result_prefetching_ set to true.");
    if (exec_params_.limit_use_default || limit_count > predicted_limit) {
        limit_count = predicted_limit;
        suppress_next_result_prefetching_ = false;
    }
    req->limit = limit_count;
}

void PgReadOp::SetRequestTotalLimit() {
    std::shared_ptr<SqlOpReadRequest> req = template_op_->request();
    // Use statement LIMIT(count + offset) for the global limit.
    int64_t limit_count = exec_params_.limit_count + exec_params_.limit_offset;
    K2LOG_D(log::pg, "Set request limit as {}", limit_count);
    req->limit = limit_count;
}

void PgReadOp::SetRowMark() {
    std::shared_ptr<SqlOpReadRequest> req = template_op_->request();

    if (exec_params_.rowmark < 0) {
        req->row_mark_type = RowMarkType::ROW_MARK_ABSENT;
    } else {
        req->row_mark_type = static_cast<RowMarkType>(exec_params_.rowmark);
    }
}

Status PgReadOp::ResetInactivePgsqlOps() {
    // Clear the existing requests.
    for (int op_index = active_op_count_; op_index < pgsql_ops_.size(); op_index++) {
        std::shared_ptr<SqlOpReadRequest> read_req = GetReadOp(op_index)->request();
        read_req->k2pgctid_column_values.clear();
        read_req->paging_state = nullptr;
    }

    // Clear row orders.
    if (batch_row_orders_.size() > 0) {
        for (int op_index = active_op_count_; op_index < pgsql_ops_.size(); op_index++) {
            batch_row_orders_[op_index].clear();
        }
    }

    return Status::OK();
}

//--------------------------------------------------------------------------------------------------

PgWriteOp::PgWriteOp(const std::shared_ptr<PgSession>& pg_session,
                        const std::shared_ptr<PgTableDesc>& table_desc,
                        const PgObjectId& table_object_id,
                        std::unique_ptr<PgWriteOpTemplate> write_op)
    : PgOp(pg_session, table_desc, table_object_id),
    write_op_(std::move(write_op)) {
}

Result<std::list<PgOpResult>> PgWriteOp::ProcessResponseImpl() {
    // Process result from doc api and check result status.
    auto result = VERIFY_RESULT(ProcessResponseResult());

    // End execution and return result.
    end_of_data_ = true;
    K2LOG_D(log::pg, "Received response for WriteOp request result size: {}", result.size());
    return result;
}

Status PgWriteOp::PopulateDmlByRowIdOps(const std::vector<std::string>& k2pgctids) {
    return STATUS(NotSupported, "PopulateDmlByRowIdOps() is not supported for PgWriteOp");
}

Status PgWriteOp::CreateRequests() {
    if (request_population_completed_) {
        return Status::OK();
    }

    // Setup a singular operator.
    // TODO: create multiple requests for multiple partitions if the information is available
    pgsql_ops_.push_back(write_op_);
    write_op_->set_active(true);
    active_op_count_ = 1;
    request_population_completed_ = true;

    return Status::OK();
}

void PgWriteOp::SetWriteTime(const uint64_t write_time) {
    write_time_ = write_time;
}

}  // namespace gate
}  // namespace k2pg

max_stars_repo_name: futurewei-cloud/chogori-sql

max_issues_repo_name: futurewei-cloud/chogori-sql

max_forks_repo_name: futurewei-cloud/chogori-sql


text: #include <boost/mpl/set/set20.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino