// ISODATA.cpp

#include "stdafx.h"
#include "ISODATA.h"
#include "Sort.h"

using namespace std;

// class ISODATA

extern int N;
extern int dim;

ISODATA::ISODATA()
{
	c = 2;	// 預期的類數；
	Nc = 1;	// 初始聚類中心個數（可以不等于c）；
	theta_n = 2;	// 每一類中允許的最少模式數目（若少于此數，就不能單獨成為一類）；
	theta_s = 1;	// 類內各分量分布的標準差上限（大于此數就分裂）；
	theta_D = 4;	// 兩類中心間的最小距離下限（若小于此數，這兩類應合并）；
	L = 1;	// 在每次迭代中可以合并的類的最多對數；
	I = 4;	// 允許的最多迭代次數；

	_d = 100;			// 總體平均距離
	Ip = 0;				// 迭代次數
	
	double D[MAXNUM][MAXNUM];	// 各類對中心間的距離

	for(int i=0;i<MAXNUM;i++)
		for(int j=0;j<MAXNUM;j++)
			D[i][j] = MAXDOUBLE;
}

ISODATA::~ISODATA()
{
}

// 設置參數
int ISODATA::SetupPattern(Pattern * pattern)
{
	for(int i=0;i<N;i++)
		x[i] = pattern[i];

	return N;
}

// 算法實現步驟
int ISODATA::Process()
{
	bool changed = true;
	// 1.預置

	// 2)將待分類的模式特征矢量x1,x2,...,xn讀入；
	// SetupPattern();

	
	// 3)選定初始聚類中心，可從待分類的模式特征矢量集{xi}中任選Nc個模式特征矢量作為初始聚類中心zj（j=1,2,...,Nc).
	InitCenter();

step1:
	// 1)設定聚類分析控制參數：
	SetupParameter();

step2:
	// 2.按最小距離原則將模式集（xi）中每個模式分到某一類中
	changed = false;
	Clustering();
	if(Ip == 0)
		cout << endl << "------------- 選取初始聚類中心 ---------------" << endl;
	else
		cout << endl << "-------------第" << Ip << "次迭代-------------" << endl;
	PrintSort();

step3:
	// 3.依據theta_n判斷合并。如果類wj中樣本數nj<theta_n，則取消該類的中心zj，Nc = Nc - 1；轉至 2.
	if(CombinBytheta_n())
		goto step2;
step4:
	// 計算分類后的參數：各類中心、類內平均距離及總體平均距離。
	CalParameter();

step5:
	// 依據Ip, Nc 判斷停止\分裂或合并。
	if(Ip == I)
	{
		theta_D = 0;
		goto step9;
	}
	else if(Nc <= c/2)
		goto step6;
	else if(Nc >= 2*c)
		goto step9;
	else if(Ip%2 == 1)	//Nc> c/2 && Nc < 2*c
		goto step6;
	else
		goto step9;

step6:
	// 計算各類類內距離的標準差矢量
	CalSigma();

step7:
	// 求出每一聚類類內距離標準差矢量sigma中的最大分量sigma_max
	// CalMaxSigma();	// 由于在Sort類中已進行計算，這里不做任何處理

step8:
	// 判斷分裂
	if(Split())
	{
		changed = true;
		Ip++;
		goto step2;
	}

step9:
	// 計算各類對中心間的距離
	CalCenterDis();

step10:
	// 依據theta_D判斷合并
	if(CombinBytheta_D())
		changed = true;

step11:
	// 判斷循環還是退出
	if(Ip >= I)
	{
		cout << endl << endl << "==========經過" << Ip << "次迭代，達到迭代次數===========" << endl;
		goto over;
	}
	else if(changed == false)
	{
		Ip++;
		cout << endl << endl << "==========經過" << Ip << "次迭代，算法收斂===========" << endl;
		goto over;
	}
	else
	{
		Ip++;

		char ch;
		cout << "本次迭代完成，是否需要改變參數（Y/N）?: ";
		cin >> ch;
		if(ch == 'y' || ch == 'Y')
			goto step1;
		else 
			goto step2;
	}

over:	
	PrintSort();
	return Ip;
}

// 3)選定初始聚類中心
void ISODATA::InitCenter()
{
	srand( (unsigned)time( NULL ) );
	int num = N;
	int no[MAXNUM];

	for(int i=0;i<N;i++)
		no[i] = i;

	for(int j=0;j<Nc;j++)
	{		
		int k = rand()%num;
		w[j].z = x[no[k]];
		w[j].z.n = 0;
		for(int l=k;l<num;l++)
			no[l] = no[l+1];
		num--;
	}
}

// step1:	設定聚類分析控制參數
void ISODATA::SetupParameter()
{
	cout << "設定聚類分析控制參數:" << endl;
	cout << "預期的類數 c:";
	cin >> c;
	cout << "初始聚類中心個數（可以不等于c）Nc:";
	cin >> Nc;
	cout << "每一類中允許的最少模式數目（若少于此數，就不能單獨成為一類）theta_n:";
	cin >> theta_n;
	cout << "類內各分量分布的標準差上限（大于此數就分裂）theta_s:";
	cin >> theta_s;
	cout << "兩類中心間的最小距離下限（若小于此數，這兩類應合并）theta_D:";
	cin >> theta_D;
	cout << "在每次迭代中可以合并的類的最多對數 L：";
	cin >> L;
	cout << "允許的最多迭代次數 I：";
	cin >> I;
}

// step2:	按最小距離原則將模式集（xi）中每個模式分到某一類中
void ISODATA::Clustering()
{
	double temp = 0.0, min = MAXDOUBLE;
	int l = 0;

	for(int j=0;j<Nc;j++)
	{
		w[j].n = 0;
		w[j].z.n = 0;
	}

	for(int i=0;i<N;i++)
	{
		min = MAXDOUBLE;
		l=0;

		for(int j=0;j<Nc;j++)
		{
			temp = Pattern::Distance(x[i], w[j].z);
			if(min > temp)
			{
				min = temp;
				l = j;
			}
		}
		w[l].Insert(x[i]);
	}
}

// step3:	依據theta_n判斷合并。如果類wj中樣本數nj<theta_n，則取消該類的中心zj，Nc = Nc - 1；轉至 step2.
bool ISODATA::CombinBytheta_n()
{
	int j=0;
	do
	{
		if(w[j].n < theta_n)
		{
			for(int k=j;k<Nc;k++)
				//w[k] = w[k+1];
				w[k].z = w[k+1].z;
			Nc--;
			// 循環中跳出！？
			return true;
		}
		j++;
	}while(j < Nc);

	return false;
}

// step4:	計算分類后的參數：各類中心、類內平均距離及總體平均距離。
void ISODATA::CalParameter()
{
	_d = 0.0;

	for(int j=0;j<Nc;j++)
	{
		w[j].CalCenter();		
		_d += w[j].n * w[j].Cal_D();
	}
	_d /= N;
}

// step6:	計算各類類內距離的標準差矢量
void ISODATA::CalSigma()
{
	for(int j=0;j<Nc;j++)
		w[j].CalSigma();
}

// step7:	求出每一聚類類內距離標準差矢量sigma中的最大分量sigma_max
void ISODATA::CalMaxSigma()
{
	// 由于在Sort類中已進行計算，這里不做任何處理
}

// step8:	判斷分裂
bool ISODATA::Split()
{
	for(int j=0;j<Nc;j++)
	{
		if( ( (w[j]._d>_d) && (w[j].n>2*(theta_n+1)) ) || (Nc <= c/2) )
		{
			int i = w[j].max;
			double sigma = w[j].sigma_max;

			for(int l=Nc;l>j;l--)
				w[l].z = w[l-1].z;

			w[j].z.x[i] -= K*sigma;
			w[j+1].z.x[i] += K*sigma;

			Nc++;
			return true;
		}
	}

	return false;
}

// step9:	計算各類對中心間的距離
void ISODATA::CalCenterDis()
{
	for(int i=0;i<Nc-1;i++)
		for(int j=i+1;j<Nc;j++)
			D[i][j] = Pattern::Distance(w[i].z, w[j].z);
}

// step10:	依據theta_D判斷合并
bool ISODATA::CombinBytheta_D()
{
	int i, j, k, l;	// 循環變量
	int num = 0;
	bool b = false;

	// 較小的類對中心
	struct
	{
		double d;
		int i;
		int j;
	}Dmin[MAXNUM];

	for(i=0;i<L;i++)
	{
		Dmin[i].d = 0.0;
		Dmin[i].i = -1;
		Dmin[i].j = -1;
	}

	// 將D[i][j]與theta_D比較，并將小于theta_D的那些D[i][j]按遞增次序排列，取前L個，Dmin[0]<Dmin[1]<Dmin[2]...
	for(i=0;i<Nc-1;i++)
		for(j=i+1;j<Nc;j++)
			if(D[i][j] < theta_D)
				for(k=0;k<L;k++)
					if(D[i][j] < Dmin[k].d)
					{
						for(l=L-1;l>k;l--)
							Dmin[l] = Dmin[l-1];
						Dmin[k].d = D[i][j];
						Dmin[k].i = i;
						Dmin[k].j = j;
						break;
					}

	// 從最小的Dmin開始，將相應的兩類合并。
	for(i=0;i<L;i++)
		if(Dmin[i].i > -1 && Dmin[i].j > -1)
			// 合并兩類
			if( w[Dmin[i].i].Combin(w[Dmin[i].j]) )
				b = true;

	// 去掉已取消的類心
	for(j=0;j<Nc;j++)
		if(w[j].z.n == -2)
		{
			for(k=j;k<Nc;k++)
				w[k].z = w[k+1].z;
			Nc--;
		}

	return b;
}

// step11:	判斷循環還是退出

// 輸出當前模式分類情況
void ISODATA::PrintSort()
{	
	cout << "共分為分為" << Nc << "類。" << endl;
	for(int i=0;i<Nc;i++)
	{
		cout << endl << "第" << i+1 << "類類心為:\t(";
		for(int j=0;j<dim-1;j++)
			cout << setw(5) << w[i].z.x[j] << ", ";
		cout << setw(5) << w[i].z.x[dim -1] << ")" << endl;

		cout << "包含的模式為:" << endl;
		for(int k=0;k<w[i].n;k++)
		{
			// Xi(12, 12, 12); 
			cout << "\tX" << w[i].x[k].n << "(";
			for(int j=0;j<dim-1;j++)
				cout << setw(5) << w[i].x[k].x[j] << ", ";
			cout << setw(5) << w[i].x[k].x[dim-1] << ");" << endl;
		}
		cout << endl;
	}
}