#include <cstdio>
#include <iostream>
#include <string>
#include <cstring>
#include <cmath>
#include <vector>
#include <queue>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include <sstream>
#include <map>
#include <random>
#include <ctime>
#include <iterator>
#include <numeric>

#include "KDTree-master/KDTree.hpp"

using namespace std;

////////////////////////////////////////////////////////////////////////////////////////// EXPERIMENT CONFIGS

const char name_str[25][205] = {
			"dataset/iris",
			"dataset/wine",
			"dataset/spam",
			"dataset/images",
			"dataset/MNIST",
			"dataset/Libras",
			"dataset/mobile",
			"dataset/zoo",
			"dataset/seeds",
			"dataset/letters",
			"dataset/phonemes"
			};
			
const char name_tag[25][205] = {
	"A",
	"B",
	"C",
	"D",
	"E",
	"F",
	"G",
	"H",
	"I",
	"J",
	"K",
};

const int num_samples[25] = {
3,
5,
793,
24,
958,
84,
112,
8,
6,
551,
396,
};
		
const double range[25][2] = {
	{0.01, 4},
	{0.01, 500},
	{0.01, 100},
	{0.01, 150},
	{0.01, 1500},
	{0.01, 1.6},
	{0.01, 800},
	{0.01, 3.5},
	{0.01, 8},
	{0.01, 10},
	{0.01, 45}
};

////////////////////////////////////////////////////////////////////////////////////////// UTILS

using Point = vector< double >;
using Points = vector< Point >;

Points points;

const double epsi = 1e-9;
const double oo = 1e100;

const int minPts = 10;

int n = 0;
int d = 0;

inline double dis(const Point& a, const Point &b)
{
	double ans = 0;
	
	for(int i = 0; i < d; i++)
		ans += (a[i] - b[i]) * (a[i] - b[i]);
		
	return sqrt(ans); 
}

inline double sq_dis(const Point& a, const Point &b)
{
	double ans = 0;
	
	for(int i = 0; i < d; i++)
		ans += (a[i] - b[i]) * (a[i] - b[i]);
		
	return ans;
}

////////////////////////////////////////////////////////////////////////////////////////// Generate Sythetic Data

const int num_centers = 3;

void generate_sythetic(int size, int dim)
{
	n = size, d = dim;
	points.clear();
	points.resize(n);
	
	clock_t time_length = clock();
	printf("*******************************************\n");
	cerr << "*******************************************\n";
	printf("Generating Sythetic Dataset\n");
	printf("size %d, dim %d\n", n, d);
	cerr << "n " << n << " d " << d << endl;
	printf("*******************************************\n");
	cerr << "*******************************************\n";
	
	default_random_engine generator;
	normal_distribution<double> distribution1(0, 10);
	normal_distribution<double> distribution2(0, 1);
	
	Points centers(num_centers);
	int seg = (int) ceil((double)n / num_centers);
	
	for(int i = 0; i < num_centers; i++)
	{
		centers[i].resize(d, 0);
		for(int j = 0; j < d; j++)
			centers[i][j] = distribution1(generator);
	}
	
	for(int i = 0; i < num_centers; i++)
	{
		for(int j = i * seg; j < n && j < (i + 1) * seg; j++)
		{
			Point x(d , 0);
			for(int k = 0; k < d; k++)
				x[k] = centers[i][k] + distribution2(generator);
			
			points[j] = x;
		}
	}
	
	time_length = clock() - time_length;
	double seconds_needed = (double)time_length / CLOCKS_PER_SEC;
	printf("*******************************************\n");
	cerr << "*******************************************\n";
	printf("Generated! Generating time %lf\n", seconds_needed);
	cerr << "Generated! Generating time " << seconds_needed << endl;
	printf("*******************************************\n");
	cerr << "*******************************************\n";
}

////////////////////////////////////////////////////////////////////////////////////////// READ DATA 

void read_data(const string &data_name)
{
	const int buffer_size = 10000000;
	static char buf[buffer_size] = {0};

	n = d = 0;
	points.clear();
	
	string points_name = data_name + "/input.txt";
	string cmd = "cp " + data_name + "/label.txt" + " label.txt";
	
	system(cmd.c_str());
	
	FILE *fp = fopen(points_name.c_str(), "r");
	
	for(;fgets(buf, buffer_size - 1, fp) != NULL; n++)
	{
		stringstream sin(buf);
		Point x;
		
		if(d == 0)
		{
			int i = 0;
			double y;
			for(; sin >> y; i++) x.push_back(y);
			d = i;
		}
		else
		{
			x.resize(d);
			for(int i = 0; sin >> x[i]; i++);
		}
		
		points.push_back(x);
	}
		
	printf("*******************************************\n");
	cerr << "*******************************************\n";
	printf("%s: %d data points loaded, dimension is %d\n", data_name.c_str(), n, d);
	cerr << "n " << n << " d " << d << endl;
	cerr << "n " << points.size() << " d " << points[0].size() << endl;
	printf("*******************************************\n");
	cerr << "*******************************************\n";
	
	fflush(stdout);
	fclose(fp);
}

////////////////////////////////////////////////////////////////////////////////////////// UNION FIND FOR CONNECTED COMPONENT 

//clustering resuts, the points in the same cluster have the same father
vector<int> father;

int find(int x)
{
	if(father[x] == x)
		return father[x];
	return father[x] = find(father[x]);
}

////////////////////////////////////////////////////////////////////////////////////////// OUTPUT AND EVALUATION  

double ARI = 0;
double AMI = 0;

void evaluate()
{
	ARI = AMI = 0;
	
	system("python3 evaluation.py > eval_tmp.txt");
	
	FILE *fp = fopen("eval_tmp.txt", "r");
	
	fscanf(fp, "%lf%lf", &ARI, &AMI);
	
	fclose(fp);
}

bool evalute_sythetic()
{
	int seg = (int) ceil((double)n / num_centers);
	
	for(int i = 1; i < n; i++)
	{
		if(i % seg == 0)
		{
			if(find(i) == find(i - 1))
				return false;
		}
		else
		{
			if(find(i) != find(i - 1))
				return false;
		}
	}
	
	for(int i = 0; i < num_centers; i++)
		for(int j = i + 1; j < num_centers; j++)
			if(find(i * seg) == find(j * seg))
				return false;

	return true;
}

void output()
{
	FILE *fo = fopen("output.txt", "w");
	
	for(int i = 0; i < n; i++)
		fprintf(fo, "%d\n", father[i]);
	
	fclose(fo);
}

////////////////////////////////////////////////////////////////////////////////////////// PARAMETER TUNING

void DBSCAN_naive(int, double, bool use_kdtree = false);
void DBSCAN(int, double, bool use_kdtree = false);
void DBSCAN_plus(int, double, int, bool, bool use_kdtree = false);
void new_DBSCAN(int, double);

void DBSCAN_tunning_linear(const string& method, const char *tag, const int itr,
		const double lb, const double ub, const int m, const int num_repeat, bool use_kdtree)
{
	double best_ARI = 0, best_ARI_T = 0, best_AMI = 0, best_AMI_T = 0, best_std_ARI = 0, best_std_AMI = 0, 
		best_ARI_time = 0, best_AMI_time = 0, best_ARI_std_time = 0, best_AMI_std_time = 0;
	int best_ARI_k = 0, best_AMI_k = 0;
	
	double T;
	int k = minPts, counter = 0;
	double add_factor = (ub - lb) / itr;
	
	cout << "###########################################\n";
	cerr << "Tuning " << method << " on dataset " << tag << " from range " << lb << " to " << ub 
		<< " with itr " << itr << " and rpt " << num_repeat <<  " use_kdtree " << use_kdtree << endl;
	cout << "Tuning " << method << " on dataset " << tag << " from range " << lb << " to " << ub 
		<< " with itr " << itr << " and rpt " << num_repeat <<  " use_kdtree " << use_kdtree << endl;
	cout << "###########################################\n";
	
	for(T = lb; T <= ub; T += add_factor)
	{
		counter++;
		
		double mean_ARI = 0, mean_second_moment_ARI = 0, std_ARI = 0;
		double mean_AMI = 0, mean_second_moment_AMI = 0, std_AMI = 0;
		double mean_time = 0, mean_second_moment_time = 0, std_time = 0;
		
		for(int t = 0; t < num_repeat; t++)
		{	
			clock_t time_length = clock();
			
			if(method == "DBSCAN")
				DBSCAN(k, T, use_kdtree);
			else if(method == "DBSCAN_naive")
				DBSCAN_naive(k, T, use_kdtree);
			else if(method == "DBSCAN_kcenter")
				DBSCAN_plus(k, T, m, true, use_kdtree);
			else if(method == "DBSCAN_uniform")	
				DBSCAN_plus(k, T, m, false, use_kdtree);
			else if(method == "new_DBSCAN")
				new_DBSCAN(k, T);
			else
			{
				printf("No method called %s\n", method.c_str());
				return;
			}
			
			time_length = clock() - time_length;
			double seconds_needed = (double)time_length / CLOCKS_PER_SEC;
			
			output();
			evaluate();
			
			
			mean_time += seconds_needed / (double)num_repeat;
			mean_ARI += ARI / (double)num_repeat;
			mean_AMI += AMI / (double)num_repeat;
			mean_second_moment_time += seconds_needed * seconds_needed / (double)num_repeat;
			mean_second_moment_ARI += ARI * ARI / (double)num_repeat;
			mean_second_moment_AMI += AMI * AMI / (double)num_repeat;
		}
		
		std_ARI = sqrt( fabs(mean_second_moment_ARI - mean_ARI * mean_ARI) );
		std_AMI = sqrt( fabs(mean_second_moment_AMI - mean_AMI * mean_AMI) );
		std_time = sqrt( fabs(mean_second_moment_time - mean_time * mean_time) );
		
		if(counter % 5 == 0)
		{
			cerr<<counter<<" times tuning."<<endl;
			cerr<<"Best ARI "<< best_ARI <<" Best AMI " << best_AMI <<" std ARI "<< best_std_ARI <<" std AMI " << best_std_AMI << endl;
			cerr<<"Best ARI time "<< best_ARI_time <<" Best AMI time " << best_AMI_time <<" best ARI std time "<< best_ARI_std_time <<" best AMI std time " << best_AMI_std_time << endl;
		}
		
		printf("%d %.4lf: ARI %.4lf, AMI %.4lf, std ARI %.4lf, std AMI %.4lf, time %lf, std time %lf\n", k, T, mean_ARI, mean_AMI, std_ARI, std_AMI, mean_time, std_time);
		fflush(stdout);
		
		if(mean_ARI > best_ARI)
			best_ARI = mean_ARI, best_ARI_k = k, best_ARI_T = T, best_std_ARI = std_ARI,
			best_ARI_time = mean_time, best_ARI_std_time = std_time;
		
		if(mean_AMI > best_AMI)
			best_AMI = mean_AMI, best_AMI_k = k, best_AMI_T = T, best_std_AMI = std_AMI,
			best_AMI_time = mean_time, best_AMI_std_time = std_time;
	}
	
	printf("===========================================\n");
	printf("Name Tag: %s, Best ARI %lf, Best AMI %lf, std ARI %lf, std AMI %lf\n", tag, best_ARI, best_AMI, best_std_ARI, best_std_AMI);
	printf("Best ARI time %lf, Best AMI time %lf, Best ARI std time %lf, Best AMI std time %lf\n", best_ARI_time, best_AMI_time, best_ARI_std_time, best_AMI_std_time);
	printf("Best ARI k %d, Best ARI eps %lf, Best AMI k %d, Best AMI eps %lf\n", best_ARI_k, best_ARI_T, best_AMI_k, best_AMI_T);
	printf("===========================================\n");
	fflush(stdout);
}


////////////////////////////////////////////////////////////////////////////////////////// ALL METHODS UNDER THIS LINE

//common members

vector< bool > is_core;
int core_points = 0;
int num_comp = 0;
vector<int> actual_core_points;

//naive DBSCAN
void DBSCAN_naive(int k, double T, bool use_kdtree)
{
	actual_core_points.clear();
	actual_core_points.resize(n, 0);
	is_core.clear();
	is_core.resize(n, false);
	father.clear();
	father.resize(n, 0);
	num_comp = 0;
	core_points = 0;
	
	for(int i = 0; i < n; i++)
	{
		father[i] = i;
		num_comp++;
	}
	
	if(use_kdtree)
	{
		KDTree tree(points);
		
		for(int i = 0; i < n; i++)
		{	
			auto ball = tree.neighborhood_indices(points[i], T);
			
			if(ball.size() >= k)
			{
				actual_core_points[core_points] = i;
				core_points++;
				is_core[i] = true;
				
				for(auto j : ball)
					father[find(j)] = find(i);
			}
		}
	}
	else
	{
		double sq_T = T * T;
		
		for(int i = 0; i < n; i++)
		{
			int ball_size = 0;
			for(int j = 0; j < n; j++)
				if( sq_dis(points[i], points[j]) <= sq_T )
					ball_size++;
			
			if(ball_size >= k)
			{
				actual_core_points[core_points] = i;
				core_points++;
				is_core[i] = true;
			}
		}
		
		for(int idx = 0 ;idx < core_points; idx++)
		{
			int i = actual_core_points[idx];
			for(int j = 0; j < n; j++)
				if(sq_dis(points[i], points[j]) <= sq_T)
					father[find(j)] = find(i);
		}	
		
	}
	
	for(int i = 0; i < n; i++)
	{
		father[i] = find(i);
		if(father[i] != i)
			num_comp--;
	}
}

//vanila DBSCAN

void DBSCAN(int k, double T, bool use_kdtree)
{
	actual_core_points.clear();
	actual_core_points.resize(n, 0);
	is_core.clear();
	is_core.resize(n, false);
	father.clear();
	father.resize(n, 0);
	num_comp = 0;
	core_points = 0;
	
	if(use_kdtree)
	{
	
		KDTree tree(points);
		
		for(int i = 0; i < n; i++)
		{
			father[i] = i;
			num_comp++;
			
			auto ball = tree.neighborhood_indices(points[i], T);
			
			if(ball.size() >= k)
			{
				actual_core_points[core_points] = i;
				core_points++;
				is_core[i] = true;
			}
		}
		
		if(core_points == 0)
			return;
		
		Points core(core_points);
		for(int idx = 0; idx < core_points; idx++)
			core[idx] = points[ actual_core_points[idx] ];
			
		KDTree core_tree(core);
		
		for(int idx_i = 0; idx_i < core_points; idx_i++)
		{
			int i = actual_core_points[idx_i];
			
			auto core_neighbor_ids = core_tree.neighborhood_indices(points[i], T);
			
			for(auto idx_j : core_neighbor_ids)
			{
				int j = actual_core_points[idx_j];
				father[find(j)] = find(i);
			}
		}
		
		for(int i = 0; i < n; i++)
			if(!is_core[i])
			{
				auto idx_j = core_tree.nearest_index(points[i]);
				int j = actual_core_points[idx_j];
				father[i] = j;
			}
	
	}
	else
	{
		double sq_T = T * T;
	
		for(int i = 0; i < n; i++)
		{
			father[i] = i;
			
			int ball_size = 0;
			for(int j = 0; j < n; j++)
				if( sq_dis(points[i], points[j]) <= sq_T )
					ball_size++;
			
			if(ball_size >= k)
			{
				actual_core_points[core_points] = i;
				core_points++;
				is_core[i] = true;
			}
		}
		
		for(int idx_i = 0; idx_i < core_points; idx_i++)
		{
			int i = actual_core_points[idx_i];
			
			for(int idx_j = idx_i + 1; idx_j < core_points; idx_j++)
			{
				int j = actual_core_points[idx_j];
				
				if( sq_dis(points[i], points[j]) <= sq_T)
						father[find(j)] = find(i);
			}
		}
		
		for(int i = 0; i < n; i++)
			if(!is_core[i])
			{
				double sq_close_core_dis = oo;
				int real_father = i;
				
				for(int idx_j = 0; idx_j < core_points; idx_j++)
				{
					int j = actual_core_points[idx_j];
					
					if(sq_dis(points[i], points[j]) <= sq_close_core_dis)
					{
						sq_close_core_dis = sq_dis(points[i], points[j]);
						real_father = j;
					}
				}
				
				father[i] = real_father;
			}
	}
	
	for(int i = 0; i < n; i++)
	{
		father[i] = find(i);
		if(father[i] != i)
			num_comp--;
	}
}

//DBSCANpp

vector< bool > is_sampled;
vector< int > per;
vector< double > min_sq_dis;

void DBSCAN_plus(int k, double T, int m, bool is_kcenter, bool use_kdtree)
{
	// all common variables
	
	actual_core_points.clear();
	actual_core_points.resize(n, 0);
	is_core.clear();
	is_core.resize(n, false);
	father.clear();
	father.resize(n, 0);
	num_comp = 0;
	core_points = 0;
	
	for(int i = 0; i < n; i++)
	{
		father[i] = i;
		num_comp++;
	}
	
	//////////////////////////////////////////////////////////////////////// 
	
	is_sampled.clear();
	is_sampled.resize(n, false);
	per.clear();
	per.resize(n, 0);
	min_sq_dis.clear();
	min_sq_dis.resize(n, oo);
	
	if(is_kcenter)
	{
		per[0] = 0;
		for(int i = 1; i < m; i++)
		{
			for(int j = 0; j < n; j++)
				min_sq_dis[j] = min(sq_dis(points[ per[i - 1] ], points[j]), min_sq_dis[j]);
		
			per[i] = distance(min_sq_dis.begin(), max_element(min_sq_dis.begin(), min_sq_dis.end()));
		}
	}
	else
	{
		for(int i = 0; i < n; i++)
			per[i] = i;
		for(int i = 0; i < n; i++)
			swap(per[i], per[rand() % (i + 1)]);
	}
	
	for(int i = 0; i < m; i++)
		is_sampled[ per[i] ] = true;
	
	if(use_kdtree)
	{
		KDTree tree(points);
		
		for(int i = 0; i < n; i++)
			if(is_sampled[i])
			{
				auto ball = tree.neighborhood_indices(points[i], T);
			
				if(ball.size() >= k)
				{
					actual_core_points[core_points] = i;
					core_points++;
					is_core[i] = true;
				}
			}
		
		if(core_points == 0)
			return;
		
		Points core(core_points);
		for(int idx = 0; idx < core_points; idx++)
			core[idx] = points[ actual_core_points[idx] ];
			
		KDTree core_tree(core);
		
		for(int idx_i = 0; idx_i < core_points; idx_i++)
		{
			int i = actual_core_points[idx_i];
			
			auto core_neighbor_ids = core_tree.neighborhood_indices(points[i], T);
			
			for(auto idx_j : core_neighbor_ids)
			{
				int j = actual_core_points[idx_j];
				father[find(j)] = find(i);
			}
		}
		
		for(int i = 0; i < n; i++)
			if(!is_core[i])
			{
				auto idx_j = core_tree.nearest_index(points[i]);
				int j = actual_core_points[idx_j];
				father[i] = j;
			}
	
	}
	else
	{
		double sq_T = T * T;
		for(int i = 0; i < n; i++)
			if(is_sampled[i])
			{
				int ball_size = 0;
				for(int j = 0; j < n; j++)
					if( sq_dis(points[i], points[j]) <= sq_T )
						ball_size++;
			
				if(ball_size >= k)
				{
					actual_core_points[core_points] = i;
					core_points++;
					is_core[i] = true;
				}
			}
		
		for(int idx_i = 0; idx_i < core_points; idx_i++)
		{
			int i = actual_core_points[idx_i];
			
			for(int idx_j = idx_i + 1; idx_j < core_points; idx_j++)
			{
				int j = actual_core_points[idx_j];
				
				if( sq_dis(points[i], points[j]) <= sq_T)
						father[find(j)] = find(i);
			}
		}
		
		for(int i = 0; i < n; i++)
			if(!is_core[i])
			{
				double sq_close_core_dis = oo;
				int real_father = i;
				
				for(int idx_j = 0; idx_j < core_points; idx_j++)
				{
					int j = actual_core_points[idx_j];
					
					if(sq_dis(points[i], points[j]) <= sq_close_core_dis)
					{
						sq_close_core_dis = sq_dis(points[i], points[j]);
						real_father = j;
					}
				}
				
				father[i] = real_father;
			}
	}
	
	for(int i = 0; i < n; i++)
	{
		father[i] = find(i);
		if(father[i] != i)
			num_comp--;
	}
}

// Our DBSCAN

vector< int > cell_hash;
const int prime = 1000000007; // an arbitrary large prime number which is used to generate a hash value;
vector< int > idx;
vector< int > tmp_cnt;
vector< int > cnt;
vector< int > tmp_core_id;
vector< int > core_id;

inline int div(const Point& a, const Point& shift, double b)
{
	int hash_value = 0;

	for(int i = 0; i < d; i++)
		hash_value = hash_value * prime + (int)floor( (a[i] + shift[i]) / b);
	
	return hash_value;
}

bool cmp(int p, int q)
{
	return cell_hash[p] < cell_hash[q];
}

inline void sort_cells(const Point& shift, double side_length)
{
	if(n <= 0) return;
	fill(tmp_cnt.begin(), tmp_cnt.end(), 1);
	fill(tmp_core_id.begin(), tmp_core_id.end(), -1);
	iota(idx.begin(), idx.end(), 0);
	
	for(int i = 0; i < n; i++)
	{
		cell_hash[i] = div(points[i], shift, side_length);
	}
	
	sort(idx.begin(), idx.end(), cmp);
	
	if(is_core[ idx[0] ])
		tmp_core_id[0] = idx[0];
	for(int i = 1; i < n; i++)
	{
		if(cell_hash[ idx[i - 1] ] == cell_hash[ idx[i] ])
			tmp_cnt[i] = tmp_cnt[i - 1] + 1, tmp_core_id[i] = tmp_core_id[i - 1];
		if(is_core[ idx[i] ])
			tmp_core_id[i] = idx[i];
	}
	
	cnt[ idx[n - 1] ] = tmp_cnt[ n - 1 ];
	core_id[ idx[n - 1] ] = tmp_core_id[ n - 1 ];
	for(int i = n - 2; i >= 0; i--)
	{
		if(cell_hash[ idx[i] ] == cell_hash[ idx[i + 1] ])
			tmp_cnt[i] = tmp_cnt[i + 1], tmp_core_id[i] = tmp_core_id[i + 1];
		cnt[ idx[i] ]  = tmp_cnt[i];
		core_id[ idx[i] ] = tmp_core_id[i];
	}
}

void new_DBSCAN(int k, double T)
{	
	// all common variables
	actual_core_points.clear();
	actual_core_points.resize(n, 0);
	is_core.clear();
	is_core.resize(n, false);
	father.clear();
	father.resize(n, 0);
	iota(father.begin(), father.end(), 0);
	num_comp = n;
	core_points = 0;
	
	//////////////////////////////////////////////////////////////////////// 
	
	int num_shifts = 2 * log(n);
	
	cell_hash.clear();
	cell_hash.resize(n, 0);
	idx.clear();
	idx.resize(n, 0);
	iota(idx.begin(), idx.end(), 0);
	cnt.clear();
	cnt.resize(n, 0);
	tmp_cnt.clear();
	tmp_cnt.resize(n, 0);
	core_id.clear();
	core_id.resize(n, 0);
	tmp_core_id.clear();
	tmp_core_id.resize(n, 0);
	
	Point shift(d, 0);
	
	//determine core
	
	sort_cells(shift, 2 * T);
	
	for(int i = 0; i < n; i++)
		if(cnt[i] >= k)
		{
			actual_core_points[core_points] = i;
			core_points++;
			is_core[ i ] = true;
		}
	
	
	//connect core points
	
	for(int rounds = 0; rounds < num_shifts; rounds++)
	{
		for(int i = 0; i < d; i++)
			shift[i] = ((double)rand()/((double)RAND_MAX + 1)) * 2 * T;
			
		sort_cells(shift, 2 * T);
		
		double sq_T = T * T;
		for(int i = 0; i < n; i++)
			if(is_core[ i ])
			{
				if(sq_dis(points[ core_id[i] ], points[ i ]) <= sq_T)
					father[find(i)] = find(core_id[i]);
			}
	}
	
	//connect non core points
	
	double side_length = T;
	bool flag_unfinished = true;
	int num_trials = num_shifts;
	for(;flag_unfinished && core_points > 0; side_length *= 2)
	{
		double sq_side_length = side_length * side_length;
		for(int rounds = 0; rounds < num_trials; rounds++)
		{
			for(int i = 0; i < d; i++)
				shift[i] = ((double)rand()/((double)RAND_MAX + 1)) * 2 * side_length;
			
			sort_cells(shift, 2 * side_length);
			
			bool side_length_too_small = true;
			
			for(int i = 0; i < n; i++)
				if(!is_core[i] && father[i] == i && core_id[i] >= 0)
				{
					if(sq_dis(points[core_id[i]], points[i]) <= sq_side_length)
					{
						father[i] = find(core_id[i]);
						side_length_too_small = false;
					}
				}
			
			if(side_length_too_small)
				break;
		}
		
		flag_unfinished = false;
		for(int i = 0; i < n; i++)
			if(!is_core[i] && father[i] == i)
			{
				flag_unfinished = true;
				break;
			}
	}
	
	// return results
	
	for(int i = 0; i < n; i++)
	{
		father[i] = find(i);
		if(father[i] != i)
			num_comp--;
	}
}

//////////////////////////////////////////////////////////////////////////////////////////

int main(int argc, char** argv)
{	
	bool run_real = true;

	string method = "new_DBSCAN";
	int num_repeats = 10;
	int num_tune = 100;
	bool use_kdtree = true;
	
	if(argc > 1)
		method = argv[1];
	if(argc > 2)
		use_kdtree = false;
	
	if(run_real)
	{
	
		vector<int> candidate_t{0,1,2,3,4,5,6,7,8,9,10};
		
		if(!(method == "DBSCAN" || method == "DBSCAN_naive" || method == "DBSCAN_kcenter"|| method == "DBSCAN_uniform" || method == "new_DBSCAN"))
			return 0;
		
		for(auto t : candidate_t)
		{
			printf("Dataset: %s\n", name_str[t]);
			
			read_data(name_str[t]);
			DBSCAN_tunning_linear(method, name_tag[t], num_tune, range[t][0], range[t][1], num_samples[t], num_repeats, use_kdtree);
		}
	
	}
	else
	{
		if(!(method == "new_DBSCAN" || method == "DBSCAN_kcenter" || method == "DBSCAN_uniform"))
			return 0;
	
		vector<int> candidate_s{100000, 200000, 300000, 400000, 500000, 600000, 700000, 800000, 10000000, 100000000};
		
		for(int s : candidate_s)
		{
		
			generate_sythetic(s, 10);
			
			int m = (int) (0.1 * pow(n, (double)d/(d+4.) ));
			clock_t time_length = clock();
			if(method == "new_DBSCAN")
				new_DBSCAN(10, 3 * sqrt(d));
			else if(method == "DBSCAN_kcenter")
				DBSCAN_plus(10, 3 * sqrt(d), m, true, use_kdtree);
			else
				DBSCAN_plus(10, 3 * sqrt(d), m, false, use_kdtree);
			cout<<"num clusters " << num_comp << endl;
			time_length = clock() - time_length;
			double seconds_needed = (double)time_length / CLOCKS_PER_SEC;
			printf("Clustering takes %lf seconds\n", seconds_needed);
			cerr<<"Clustering takes  "<< seconds_needed << " seconds\n";
			cerr<<"Start evaluation\n";
			if(evalute_sythetic())
			{
				printf("Pefect clustering!\n=====================================\n");
			}
			else
			{
					printf("Fail to get perfect clustering!\n=====================================\n");
					cerr<<"Fail to get perfect clustering!\n=====================================\n";
			}
			fflush(stdout);
		}
		
		vector<int> candidate_D{10, 100, 300};
		int s = 1000000;
		for(int D : candidate_D)
		{
			generate_sythetic(s, D);
			
			int m = (int) (0.1 * pow(n, (double)d/(d+4.) ));
			clock_t time_length = clock();
			if(method == "new_DBSCAN")
				new_DBSCAN(10, 3 * sqrt(d));
			else if(method == "DBSCAN_kcenter")
				DBSCAN_plus(10, 3 * sqrt(d), m, true, use_kdtree);
			else
				DBSCAN_plus(10, 3 * sqrt(d), m, false, use_kdtree);
			cout<<"num clusters " << num_comp << endl;
			time_length = clock() - time_length;
			double seconds_needed = (double)time_length / CLOCKS_PER_SEC;
			printf("Clustering takes %lf seconds\n", seconds_needed);
			cerr<<"Clustering takes  "<< seconds_needed << " seconds\n";
			cerr<<"Start evaluation\n";
			if(evalute_sythetic())
			{
				printf("Pefect clustering!\n=====================================\n");
			}
			else
			{
					printf("Fail to get perfect clustering!\n=====================================\n");
					cerr<<"Fail to get perfect clustering!\n=====================================\n";
			}
			fflush(stdout);
		}
	}
	
	return 0;
}




















