//Coreset Construction

#include "coreset.h"


parti group(const parti& _Light, const datapoint& c, const double& err, int z) {
	parti _G;
	dataset C;
	double cur = 0;
	for (dataset P : _Light) {
		if (cur + cost(P, c, z) <= err) {
			Union(C, P);
			cur += cost(P, c, z);
		}
		else {
			_G.push_back(C);
			C = P;
			cur = cost(P, c, z);
		}
	}
	_G.push_back(C);
	return _G;
}


pair<parti, parti> WZ_Decomposition(dataset X, const datapoint& c, const double& err, int z) {
	parti _Rings, _W, _Z;
	double maxd = 1, mind = 1e18;
	for (datapoint x : X) {
		maxd = max(maxd, dist(x, c));
		if (dist(x, c) > 0) mind = min(mind, dist(x, c));
	}
	_Rings.resize(int(log2(maxd / mind)) + 2);
	for (datapoint x : X) {
		if (dist(x, c) == 0) _Rings[0].push_back(x);
		else _Rings[int(log2(dist(x, c) / mind)) + 1].push_back(x);
	}
	parti _Light;
	for (dataset P : _Rings) {
		if (cost(P, c, z) > err) {
			_W.push_back(P);
			Union(_Z, group(_Light, c, err, z));
			_Light.clear();
		}
		else {
			_Light.push_back(P);
		}
	}
	Union(_Z, group(_Light, c, err, z));
	return make_pair(_W, _Z);
}

//X: group
dataset Two_point_Coreset(dataset X, const datapoint& c) {
	if (X.empty()) return X;
	datapoint p_c = X[0], p_f = X[0];
	for (datapoint x : X) {
		if (dist(x, c) < dist(p_c, c)) p_c = x;
		if (dist(x, c) > dist(p_f, c)) p_f = x;
	}
	double d_c = dist(p_c, c), d_f = dist(p_f, c);
	if (d_c == d_f) {
		p_c[0] = 0;
		for (datapoint x : X) p_c[0] += x[0];
		dataset temp; temp.push_back(p_c);
		return temp;
	}
	p_c[0] = 0;
	p_f[0] = 0;
	for (datapoint x : X) {
		p_f[0] += x[0] * (dist(x, c) - d_c) / (d_f - d_c);
		p_c[0] += x[0] * (d_f - dist(x, c)) / (d_f - d_c);
	}
	dataset temp; temp.push_back(p_c); temp.push_back(p_f);
	return temp;
}


dataset Uniform_Coreset(dataset X, const int& n) {
	dataset Coreset;
	vector<double> distri;
	double sum = 0;
	for (int i = 0; i < X.size(); i++) {
		distri.push_back(X[i][0]);
		sum += X[i][0];
	}
	sampler sp; sp.init(distri);
	for (int i = 1; i <= n; i++) {
		datapoint x = X[sp.sample()];
		x[0] = sum / n;
		Coreset.push_back(x);
	}
	return Coreset;
}

dataset HJLW23_Coreset(dataset X, dataset C_approx, double m, int x, int z, double thr) {
	dataset Coreset;
	double err;
	vector<pair<double, int>> d;//(dist,index)//ascending order
	for (int i = 0; i < X.size(); i++) {
		d.push_back({ dist(X[i], C_approx), i });
	}
	sort(d.begin(), d.end());
	dataset Y = X;
	for (int i = 0; i < X.size(); i++) X[i] = Y[d[i].second];
	for (int i = X.size(); ~i; i--) {
		if (m < 1e-7) {
			X.resize(i); break;
		}
		m -= X[i - 1][0];
		Coreset.push_back(X[i - 1]);
	}
	if (x == 0) return Coreset;
	parti _Clusters;
	_Clusters.resize(C_approx.size());
	for (int i = 0; i < (int)C_approx.size(); i++) C_approx[i][0] = i;
	for (datapoint x : X) {
		_Clusters[NN(x, C_approx)[0]].push_back(x);
	}
	parti _H;
	for (int i = 0; i < (int)C_approx.size(); i++) {
		if (_Clusters[i].empty()) continue;
		auto D = WZ_Decomposition(_Clusters[i], C_approx[i], thr * cost(_Clusters[i], C_approx[i], z), z);
		parti _W = D.first, _Z = D.second;
		for (dataset S : _Z) {
			dataset tpCoreS = Two_point_Coreset(S, C_approx[i]);
			Union(Coreset, tpCoreS);
			x -= tpCoreS.size();
		}
		Union(_H, _W);
	}
	if (x < 0) {
		cerr << "Coreset construction failed." << endl;
		return Coreset;
	}
	if (_H.empty()) return Coreset;
	
	vector<int> sz;
	for (int i = 0; i < _H.size(); i++) sz.push_back(x / (int)_H.size());
	sampler sp; sp.init_uniform(_H.size());
	for (int i = 0; i < x % (int)_H.size(); i++) sz[sp.sample()]++;
	for (int i = 0; i < _H.size(); i++) {
		Union(Coreset, Uniform_Coreset(_H[i], sz[i]));
	}
	return Coreset;
}


dataset HLLW25_Coreset(dataset X, dataset C_approx, double m, int x, int z, double thr) {
	dataset Coreset;
	double err;
	vector<pair<double, int>> d;//(dist,index)//ascending order
	for (int i = 0; i < X.size(); i++) {
		d.push_back({ dist(X[i], C_approx), i });
	}
	sort(d.begin(), d.end());
	dataset Y = X;
	for (int i = 0; i < X.size(); i++) X[i] = Y[d[i].second];
	for (int i = X.size(); ~i; i--) {
		if (m < 1e-7) {
			X.resize(i); break;
		}
		m -= X[i - 1][0];
		Coreset.push_back(X[i - 1]);
	}
	if (x == 0) return Coreset;
	parti _Clusters;
	_Clusters.resize(C_approx.size());
	for (int i = 0; i < (int)C_approx.size(); i++) C_approx[i][0] = i;
	for (datapoint x : X) {
		_Clusters[NN(x, C_approx)[0]].push_back(x);
	}


	vector<parti> _H;
	for (int i = 0; i < (int)C_approx.size(); i++) {// i th cluster
		if (_Clusters[i].empty()) continue;
		auto D = WZ_Decomposition(_Clusters[i], C_approx[i], thr * cost(_Clusters[i], C_approx[i], z), z);
		parti _W = D.first, _Z = D.second;
		for (dataset S : _Z) {
			dataset tpCoreS = Two_point_Coreset(S, C_approx[i]);
			Union(Coreset, tpCoreS);
			x -= tpCoreS.size();
		}
		_H.push_back(_W);
	}


	if (x < 0) {
		cerr << "Coreset construction failed." << endl;
		return Coreset;
	}
	if (_H.empty()) return Coreset;

	int cx = x / _H.size();
	sampler sp1; sp1.init_uniform(_H.size());
	vector<int> cluster_sz;
	for (int i = 0; i < _H.size(); i++) { cluster_sz.push_back(cx); }
	for (int i = 0; i < x % _H.size(); i++) cluster_sz[sp1.sample()]++;

	for (int i = 0; i < _H.size(); i++)
	{
		parti _W = _H[i];// i th ring set
		double clusterCost = cost(_Clusters[i], C_approx[i], z);
		int remainSize = cluster_sz[i];
		vector<int> sz;
		for (int j = 0; j < _W.size(); j++)// each ring
		{
			double ringCost = cost(_W[j], C_approx[i], z);
			double lambdaR = ringCost / clusterCost;
			int sampleSize = cluster_sz[i] * lambdaR;
			sz.push_back(sampleSize);
			remainSize -= sampleSize;
		}
		sampler sp; sp.init_uniform(_W.size());
		for (int j = 0; j < remainSize; j++) sz[sp.sample()]++;
		for (int j = 0; j < _W.size(); j++) {
			Union(Coreset, Uniform_Coreset(_W[j], sz[j]));
		}
	}
	return Coreset;
}

datapoint MeanPoint(dataset bucket)
{
	double sum = 0;
	for (int i = 0; i < bucket.size(); i++)
	{
		sum += bucket[i][1];
	}
	sum = sum / bucket.size();
	datapoint x; x.resize(2);
	x[0] = 1; x[1] = sum;
	return x;
}

//Our high dimension method
dataset Our_Coreset(dataset X, dataset C_approx, double m, int x, int z, double thr) {
	dataset Coreset;
	dataset Lstar;
	double err;
	vector<pair<double, int>> d;//(dist,index)//ascending order
	for (int i = 0; i < X.size(); i++) {
		d.push_back({ dist(X[i], C_approx), i });
	}
	sort(d.begin(), d.end());
	dataset Y = X;
	for (int i = 0; i < X.size(); i++) X[i] = Y[d[i].second];
	for (int i = X.size(); ~i; i--) {
		if (m < 1e-7) {
			X.resize(i); break;
		}
		m -= X[i - 1][0];
		Lstar.push_back(X[i - 1]);
	}
	if (x == 0) return Coreset;
	parti _Clusters;
	_Clusters.resize(C_approx.size());
	for (int i = 0; i < (int)C_approx.size(); i++) C_approx[i][0] = i;
	for (datapoint x : X) {
		_Clusters[NN(x, C_approx)[0]].push_back(x);
	}


	vector<parti> _H;
	for (int i = 0; i < (int)C_approx.size(); i++) {// i th cluster
		if (_Clusters[i].empty()) continue;
		auto D = WZ_Decomposition(_Clusters[i], C_approx[i], thr * cost(_Clusters[i], C_approx[i], z), z);
		parti _W = D.first, _Z = D.second;
		for (dataset S : _Z) {
			dataset tpCoreS = Two_point_Coreset(S, C_approx[i]);
			Union(Coreset, tpCoreS);
			x -= tpCoreS.size();
		}
		_H.push_back(_W);
	}


	if (x < 0) {
		cerr << "Coreset construction failed." << endl;
		return Coreset;
	}
	if (_H.empty()) return Coreset;

	int LstarSize = x / (_H.size() + 1);
	x = x - LstarSize;
	int cx = x / _H.size();
	sampler sp1; sp1.init_uniform(_H.size());
	vector<int> cluster_sz;
	for (int i = 0; i < _H.size(); i++) { cluster_sz.push_back(cx); }
	for (int i = 0; i < x % _H.size(); i++) cluster_sz[sp1.sample()]++;

	for (int i = 0; i < _H.size(); i++)
	{
		parti _W = _H[i];// ring set of the i th cluster
		double clusterCost = cost(_Clusters[i], C_approx[i], z);
		int remainSize = cluster_sz[i];
		vector<int> sz;
		for (int j = 0; j < _W.size(); j++)// each ring
		{
			double ringCost = cost(_W[j], C_approx[i], z);
			double lambdaR = ringCost / clusterCost;
			int sampleSize = cluster_sz[i] * lambdaR;
			sz.push_back(sampleSize);
			remainSize -= sampleSize;
		}
		sampler sp; sp.init_uniform(_W.size());
		for (int j = 0; j < remainSize; j++) sz[sp.sample()]++;
		for (int j = 0; j < _W.size(); j++) {
			Union(Coreset, Uniform_Coreset(_W[j], sz[j]));
		}
	}
	Union(Coreset, Uniform_Coreset(Lstar, LstarSize));
	return Coreset;
}

dataset Coreset1d(dataset X, double eps)
{
	int n = X.size();

	double OPT = cost(X, X[floor(n / 2)], 1);
	vector<pair<double, int>> d;//(dist,index)//ascending order
	for (int i = 0; i < X.size(); i++) {
		d.push_back({ X[i][1], i });
	}
	sort(d.begin(), d.end());
	int x = floor(eps * n);
	dataset B_L; dataset B_MR; dataset B_ML; dataset B_R;
	dataset coreset;
	for (int i = 0; i < d.size(); i++)
	{
		if (i <= x - 1) B_L.push_back(X[d[i].second]);
		else
		{
			if (i <= n / 2) B_ML.push_back(X[d[i].second]);
			else
				if (i <= n - x - 1) B_MR.push_back(X[d[i].second]);
				else B_R.push_back(X[d[i].second]);
		}
	}

	if (!B_L.empty())
	{
		datapoint SL = MeanPoint(B_L); SL[0] = B_L.size(); coreset.push_back(SL);
	}
	if (!B_R.empty())
	{
		datapoint SR = MeanPoint(B_R); SR[0] = B_R.size(); coreset.push_back(SR);
	}


	//B_ML processing
	vector<double> prefix1;

	double cost1 = cost(X, B_ML[0], 1);
	int leftBML = x + 1;
	int curRight = X.size() - leftBML;
	prefix1.push_back(cost1); // Total cost at each point
	for (int i = 1; i < B_ML.size(); i++)
	{
		cost1 = cost1 + (leftBML - (curRight)) * (B_ML[i][1] - B_ML[i - 1][1]);
		prefix1.push_back(cost1);
		curRight -= 1;
		leftBML += 1;
	}
	vector<pair<int, int> > costLs;//(block num,index in B_L)
	int n1 = B_ML.size();
	for (int i = 0; i < n1; i++)
	{
		double curCost = prefix1[i];
		int num = floor(log2(curCost / OPT));
		costLs.push_back({ num,i });
	}
	sort(costLs.begin(), costLs.end());


	int curBlock = costLs[0].first;
	int curBlockIndex = 0;
	double uppercost = eps * pow(2, costLs[0].first) * OPT;

	//Greedy Part: partition block into buckets
	double preMu = B_ML[costLs[0].second][1];
	int nBlock = 1; // muber of points in current bucket
	int indexAfterMu = B_ML.size();
	for (int i = 1; i < B_ML.size(); i++)
	{
		if (B_ML[i][1] > preMu) {
			indexAfterMu = i;
			break;
		}
	}

	double preCu = 0;

	int i = 1;
	if (costLs.size() == 1)
	{
		datapoint mu = B_ML[0]; mu[0] = 1;
		coreset.push_back(mu);
	}

	while (i < costLs.size())
	{
		dataset bucket;
		while (costLs[i].first == curBlock)
		{
			double curValue = B_ML[costLs[i].second][1];
			bucket.push_back(B_ML[costLs[i].second]);
			double curMu = (preMu * nBlock + curValue) / (nBlock + 1);
			int indexAfterCurMu = indexAfterMu;
			double sumD = 0;

			if (indexAfterMu < B_ML.size() && curMu != preMu)
			{
				indexAfterCurMu = B_ML.size();
				for (int j = indexAfterMu; j < B_ML.size(); j++)
				{
					if (B_ML[j][1] > curMu) {
						indexAfterCurMu = j;
						break;
					}
					else
					{
						sumD = sumD + (curMu - B_ML[j][1]) - (B_ML[j][1] - preMu);
					}

				}
			}
			// when i is added to the bucket
			double cumulative = preCu + (min(indexAfterMu - curBlockIndex, i - curBlockIndex + 1) - max(0, i - indexAfterCurMu)) * (curMu - preMu) + sumD + curValue - curMu;

			if (cumulative >= uppercost)
			{
				bucket.erase(bucket.end() - 1);
				datapoint mu; mu.resize(2); mu[0] = bucket.size(); mu[1] = preMu;
				coreset.push_back(mu);
				bucket.clear();
				bucket.push_back(B_ML[costLs[i].second]);
				curBlockIndex = i;
				preMu = B_ML[costLs[i].second][1];
				indexAfterMu = B_ML.size();
				for (int j = i; j < B_ML.size(); j++)
				{
					if (B_ML[j][1] > preMu) {
						indexAfterMu = j;
						break;
					}
				}
				nBlock = 1;
				preCu = 0;
			}
			else
			{
				indexAfterMu = indexAfterCurMu;
				preMu = curMu;
				preCu = cumulative;
				nBlock++;
			}

			i++;
			if (i >= costLs.size())
			{
				datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
				coreset.push_back(mu);
				break;
			}

		}

		//Get into a new block
		if (i < costLs.size())
		{
			datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
			coreset.push_back(mu);
			curBlock = costLs[i].first;
			uppercost = eps * pow(2, costLs[i].first) * OPT;
			curBlockIndex = i;
			preMu = B_ML[costLs[i].second][1];
			nBlock = 1;
			preCu = 0;
			indexAfterMu = B_ML.size();
			for (int j = i; j < B_ML.size(); j++)
			{
				if (B_ML[j][1] > preMu) {
					indexAfterMu = j;
					break;
				}
			}
		}

	}

	//B_MR processing
	prefix1.clear();

	cost1 = cost(X, B_MR[0], 1);
	int leftBMR = x + 1 + B_MR.size();
	curRight = X.size() - leftBMR;
	prefix1.push_back(cost1);
	for (int i = 1; i < B_MR.size(); i++)
	{
		cost1 = cost1 + (leftBMR - (curRight)) * (B_MR[i][1] - B_MR[i - 1][1]);
		prefix1.push_back(cost1);
		curRight -= 1;
		leftBMR += 1;
	}
	vector<pair<int, int> > costRs;//(block num,index in B_L)
	n1 = B_MR.size();
	for (int i = 0; i < n1; i++)
	{
		double curCost = prefix1[i];
		int num = floor(log2(curCost / OPT));
		costRs.push_back({ num,i });
	}
	sort(costRs.begin(), costRs.end());


	curBlock = costRs[0].first;
	curBlockIndex = 0;
	uppercost = eps * pow(2, costRs[0].first) * OPT;

	//Greedy Part: partition block into buckets
	preMu = B_MR[costRs[0].second][1];
	nBlock = 1;
	indexAfterMu = B_MR.size();
	for (int i = 1; i < B_MR.size(); i++)
	{
		if (B_MR[i][1] > preMu) {
			indexAfterMu = i;
			break;
		}
	}

	preCu = 0;
	i = 1;
	if (costRs.size() == 1)
	{
		datapoint mu = B_MR[0]; mu[0] = 1;
		coreset.push_back(mu);
	}

	while (i < costRs.size())
	{
		dataset bucket;
		while (costRs[i].first == curBlock)
		{
			double curValue = B_MR[costRs[i].second][1];
			bucket.push_back(B_MR[costRs[i].second]);
			double curMu = (preMu * nBlock + curValue) / (nBlock + 1);
			int indexAfterCurMu = indexAfterMu;
			double sumD = 0;

			if (indexAfterMu < B_MR.size() && curMu != preMu)
			{
				indexAfterCurMu = B_MR.size();
				for (int j = indexAfterMu; j < B_MR.size(); j++)
				{
					if (B_MR[j][1] > curMu) {
						indexAfterCurMu = j;
						break;
					}
					else
					{
						sumD = sumD + (curMu - B_MR[j][1]) - (B_MR[j][1] - preMu);
					}

				}
			}
			double cumulative = preCu + (min(indexAfterMu - curBlockIndex, i - curBlockIndex + 1) - max(0, i - indexAfterCurMu)) * (curMu - preMu) + sumD + curValue - curMu;

			if (cumulative >= uppercost)
			{
				bucket.erase(bucket.end() - 1);
				datapoint mu; mu.resize(2); mu[0] = bucket.size(); mu[1] = preMu;
				coreset.push_back(mu);
				bucket.clear();
				bucket.push_back(B_MR[costRs[i].second]);
				curBlockIndex = i;
				preMu = curValue;
				nBlock = 1;
				preCu = 0;
				indexAfterMu = B_MR.size();
				for (int j = i; j < B_MR.size(); j++)
				{

					if (B_MR[j][1] > preMu) {
						indexAfterMu = j;
						break;
					}
				}

			}
			else
			{
				indexAfterMu = indexAfterCurMu;
				preMu = curMu;
				preCu = cumulative;
				nBlock++;
			}

			i++;
			if (i >= costRs.size())
			{
				datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
				coreset.push_back(mu);
				break;
			}

		}
		if (i < costRs.size())
		{
			datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
			coreset.push_back(mu);
			curBlock = costRs[i].first;
			uppercost = eps * pow(2, costRs[i].first) * OPT;
			curBlockIndex = i;
			preMu = B_MR[costRs[i].second][1];
			nBlock = 1;
			preCu = 0;
			indexAfterMu = B_MR.size();
			for (int j = i; j < B_MR.size(); j++)
			{
				if (B_MR[j][1] > preMu) {
					indexAfterMu = j;
					break;
				}
			}
		}

	}



	return coreset;
}

//Our 1 dimension 1-median method
dataset Bucket_Coreset(dataset X, double m, double eps)
{
	int n = X.size();
	vector<pair<double, int>> d;//(dist,index)//ascending order
	for (int i = 0; i < X.size(); i++) {
		d.push_back({ X[i][1], i });
	}
	sort(d.begin(), d.end());

	dataset P_L; dataset P_R; dataset P_M;
	for (int i = 0; i < d.size(); i++)
	{
		if (i <= m - 1) P_L.push_back(X[d[i].second]);
		else
		{
			if (i <= n - m - 1) P_M.push_back(X[d[i].second]);
			else P_R.push_back(X[d[i].second]);
		}
	}
	dataset coreset = Coreset1d(P_M, eps / (double)3);


	// compute the approximation c^\star
	dataset Cstar; Cstar.push_back(X[d[n / 2].second]);

	d.clear();
	for (int i = 0; i < X.size(); i++) {
		d.push_back({ dist(X[i], Cstar), i });
	}
	sort(d.begin(), d.end());
	int rmax = d[X.size() - m - 1].first;
	datapoint CL; datapoint CR; bool cl = false; bool cr = false;
	for (int i = d.size() - m - 1; i >= 0; i--)
	{
		if (X[d[i].second][1] <= Cstar[0][1] && cl == false)
		{
			CL = X[d[i].second]; cl = true;
		}
		if (X[d[i].second][1] > Cstar[0][1] && cr == false)
		{
			CR = X[d[i].second]; cr = true;
		}
		if (cl && cr) { break; }
	}
	vector<dataset> _B; dataset PLL; dataset PLR;
	for (int i = 0; i < P_L.size(); i++)
	{
		if (P_L[i][1] < CL[1]) { PLL.push_back(P_L[i]); }
		else PLR.push_back(P_L[i]);
	}
	dataset PRL; dataset PRR;
	for (int i = 0; i < P_R.size(); i++)
	{
		if (P_R[i][1] <= CR[1]) { PRL.push_back(P_R[i]); }
		else PRR.push_back(P_R[i]);
	}
	_B.push_back(PLL); _B.push_back(PLR); _B.push_back(PRL); _B.push_back(PRR);
	for (int i = 0; i < _B.size(); i++)
	{
		dataset curData = _B[i];
		if (curData.empty())
		{
			continue;
		}
		datapoint pc;
		if (i <= 1) { pc = CL; }
		else { pc = CR; }

		// d: distance to c_L(c_R)
		d.clear();
		for (int j = 0; j < curData.size(); j++) {//j's range only in curData
			d.push_back({ dist(curData[j], pc), j });
		}
		sort(d.begin(), d.end()); //dist to cstar, index in curData
		dataset Bfar;
		int nIndex = -1;
		for (int j = d.size() - 1; j >= 0; j--)
		{
			if (d[j].first >= rmax) { Bfar.push_back(curData[d[j].second]); }
			else
			{
				nIndex = j; break;
			}
		}
		int step = eps * (double)n / (double)16;
		dataset bucket;
		if (step == 0)
		{
			step = 1;
		}

		//Bfar bucket partition
		for (int j = 0; j < Bfar.size(); j++)
		{
			if (j % step == 0)
			{
				if (!bucket.empty())
				{
					datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
					coreset.push_back(mu);
				}
				bucket.clear();
			}
			bucket.push_back(Bfar[j]);
		}
		if (!bucket.empty())
		{
			datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
			coreset.push_back(mu);
		}
		if (nIndex == -1)// All points are in far buckets
		{
			continue;
		}

		// Partiton points not in B_far into blocks
		int lIndex = nIndex;
		for (int j = 0; j <= nIndex; j++) // <= nIndex points in B_0 or B_i
		{
			if (d[j].first >= eps * rmax * 2) { lIndex = j; break; }
		}

		int curBlock;
		double uppercumu;
		int j = 0;
		double preMu = 0;
		int nBlock = 1;
		int indexAfterMu;
		double preCu = 0;
		int curBlockIndex = 0;
		bool B0 = true;
		bucket.clear();
		while (j <= nIndex)
		{
			int pBlock = 0;
			if (j >= lIndex) // The points out of block 0
			{
				pBlock = log2(d[j].first / (eps * rmax));
			}
			if (j == 0)// Initialization
			{
				preMu = d[j].first;
				curBlock = pBlock;
				uppercumu = eps * eps * pow(2, curBlock) * n * rmax / (double)288;
				nBlock = 1;
				indexAfterMu = nIndex + 1;
				bucket.push_back(curData[d[j].second]);
				for (int g = 1; g <= nIndex; g++)
				{
					if (d[g].first > preMu) {
						indexAfterMu = g;
						break;
					}
				}
				j++;
				if (j > nIndex)
				{
					datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
					coreset.push_back(mu);
					break;
				}
				continue;


			}

			while (pBlock == curBlock)
			{

				double curValue = d[j].first;
				bucket.push_back(curData[d[j].second]);
				double curMu = (preMu * nBlock + curValue) / (nBlock + 1);
				int indexAfterCurMu = indexAfterMu;
				double sumD = 0;

				if (indexAfterMu <= nIndex && curMu != preMu)
				{
					indexAfterCurMu = nIndex + 1;
					for (int j = indexAfterMu; j <= nIndex; j++)
					{
						if (d[j].first > curMu) {
							indexAfterCurMu = j;
							break;
						}
						else
						{
							sumD = sumD + (curMu - d[j].first) - (d[j].first - preMu);
						}

					}
				}
				double cumulative = preCu + (min(indexAfterMu - curBlockIndex, j - curBlockIndex + 1) - max(0, j - indexAfterCurMu)) * (curMu - preMu) + sumD + curValue - curMu;

				if (cumulative > uppercumu || bucket.size() > step)
				{
					bucket.erase(bucket.end() - 1);
					datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
					coreset.push_back(mu);
					bucket.clear();
					bucket.push_back(curData[d[j].second]);
					curBlockIndex = j;
					preMu = curValue;
					nBlock = 1;
					preCu = 0;
					indexAfterMu = nIndex + 1;
					for (int g = j; g <= nIndex; g++)
					{
						if (d[g].first > preMu) {
							indexAfterMu = g;
							break;
						}
					}

				}
				else
				{
					indexAfterMu = indexAfterCurMu;
					preMu = curMu;
					preCu = cumulative;
					nBlock++;
				}

				j++;


				if (j > nIndex)
				{
					datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
					coreset.push_back(mu);
					break;
				}
				if (j >= lIndex)
				{
					pBlock = log2(d[j].first / (eps * rmax));
				}
				else
				{
					pBlock = 0;
				}

			}

			//current point not in the previous block
			if (j <= nIndex)
			{
				datapoint mu = MeanPoint(bucket); mu[0] = bucket.size();
				coreset.push_back(mu);
				curBlock = pBlock;
				uppercumu = eps * eps * pow(2, curBlock) * n * rmax / (double)288;
				curBlockIndex = j;
				preMu = d[j].first;
				nBlock = 1;
				preCu = 0;
				indexAfterMu = nIndex + 1;
				bucket.clear();
				bucket.push_back(curData[d[j].second]);
				j++;
				for (int g = j; g <= nIndex; g++)
				{
					if (d[g].first > preMu) {
						indexAfterMu = g;
						break;
					}
				}

			}
		}


	}


	return coreset;
}

dataset Outlier_Uniform(dataset X, dataset C_approx, double m, int x) {
	dataset Coreset;
	vector<pair<double, int>> d;
	for (int i = 0; i < X.size(); i++) {
		d.push_back({ dist(X[i], C_approx), i });
	}
	sort(d.begin(), d.end());
	dataset Y = X;
	for (int i = 0; i < X.size(); i++) X[i] = Y[d[i].second];
	for (int i = X.size(); i; i--) {
		if (m < 1e-7) {
			X.resize(i); break;
		}
		m -= X[i - 1][0];
		Coreset.push_back(X[i - 1]);
	}
	Union(Coreset, Uniform_Coreset(X, x));
	return Coreset;
}
