#include "utility.h"



double postCorrection_1(int* localGenoDistr, bool flag)
{
	float mu[36];
	float mu0[36] = { 0 };
	float mu_ik[6] = { 0 };
	float mu_jk[6] = { 0 };
	float mu_kt[4] = { 0 };
	float mu_ijt[18] = { 0 };
	float mu_ijk[18] = { 0 };

	float n_ik[6] = { 0 };
	float n_jk[6] = { 0 };
	float n_kt[4] = { 0 };
	float n_ijk[18] = { 0 };
	float n_ijt[18] = { 0 };

	float muError = 36.0;
	float Likelihood = 0.0;
	float tmp = 0;
	for (int i = 0; i < 36; i++)mu[i] = 1;

	for (int i = 0; i < 3; i++)
	{
		for (int k = 0; k < 2; k++)
		{
			for (int j = 0; j < 3; j++)
			{
				for (int t = 0; t < 2; t++)
				{
					n_ik[i * 2 + k] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
					n_jk[j * 2 + k] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
					n_kt[k * 2 + t] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
					n_ijt[i * 6 + j * 2 + t] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
					n_ijk[i * 6 + j * 2 + k] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
				}
			}
		}
	}
	while (muError>1)
	{
		for (int i = 0; i<36; i++)
		{
			mu0[i] = mu[i];
		}
		//step2: mu_ik and n_ik
		if (!flag)
		{
			for (int i = 0; i < 6; i++)
			{
				mu_ik[i] = 0;
				//n_ik[i] = 0;
			}
			for (int i = 0; i < 3; i++)
			{
				for (int k = 0; k < 2; k++)
				{
					for (int j = 0; j < 3; j++)
					{
						for (int t = 0; t < 2; t++)
						{
							mu_ik[i * 2 + k] += mu[i * 12 + j * 4 + k * 2 + t];
							//n_ik[i * 2 + k] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
						}
					}
				}
			}

			//mu_ijkt=mu_ijkt*n_ik/mu_ik

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t = 0; t < 2; t++)
						{
							if (mu_ik[i * 2 + k] > 0)
								mu[i * 12 + j * 4 + k * 2 + t] = mu[i * 12 + j * 4 + k * 2 + t] * n_ik[i * 2 + k] / mu_ik[i * 2 + k];
							else
								mu[i * 12 + j * 4 + k * 2 + t] = 0;
						}
					}
				}
			}

			//step3:mu_jk and n_jk
			for (int i = 0; i < 6; i++)
			{
				mu_jk[i] = 0;
				//n_jk[i] = 0;
			}
			for (int j = 0; j < 3; j++)
			{
				for (int k = 0; k < 2; k++)
				{
					for (int i = 0; i < 3; i++)
					{
						for (int t = 0; t < 2; t++)
						{
							mu_jk[j * 2 + k] += mu[i * 12 + j * 4 + k * 2 + t];
							//n_jk[j * 2 + k] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
						}
					}
				}
			}

			//mu_ijkt = mu_ijkt*n_jk/mu_jk
			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t = 0; t < 2; t++)
						{
							if (mu_jk[j * 2 + k] > 0)
								mu[i * 12 + j * 4 + k * 2 + t] = mu[i * 12 + j * 4 + k * 2 + t] * n_jk[j * 2 + k] / mu_jk[j * 2 + k];
							else
								mu[i * 12 + j * 4 + k * 2 + t] = 0;
						}
					}
				}
			}
		}
		//step4: mu_kt and n_kt
		for (int i = 0; i<4; i++)
		{
			mu_kt[i] = 0;
			//n_kt[i] = 0;
		}
		for (int k = 0; k<2; k++)
		{
			for (int t = 0; t<2; t++)
			{
				for (int i = 0; i<3; i++)
				{
					for (int j = 0; j<3; j++)
					{
						mu_kt[k * 2 + t] += mu[i * 12 + j * 4 + k * 2 + t];
						//n_kt[k * 2 + t] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
					}
				}
			}
		}

		//mu_ijkt = mu_ijkt*n_kt/mu_kt
		for (int i = 0; i<3; i++)
		{
			for (int j = 0; j<3; j++)
			{
				for (int k = 0; k<2; k++)
				{
					for (int t = 0; t<2; t++)
					{
						if (mu_kt[k * 2 + t]>0)
							mu[i * 12 + j * 4 + k * 2 + t] = mu[i * 12 + j * 4 + k * 2 + t] * n_kt[k * 2 + t] / mu_kt[k * 2 + t];
						else
							mu[i * 12 + j * 4 + k * 2 + t] = 0;
					}
				}
			}
		}
		//mu_ijkt = mu_ijkt*n_ijt/mu_ijt
		for (int i = 0; i<18; i++)
		{
			mu_ijt[i] = 0;
			//n_ijt[i] = 0;
		}
		for (int i = 0; i<3; i++)
		{
			for (int j = 0; j<3; j++)
			{
				for (int t = 0; t<2; t++)
				{
					for (int k = 0; k<2; k++)
					{
						mu_ijt[i * 6 + j * 2 + t] += mu[i * 12 + j * 4 + k * 2 + t];
						//n_ijt[i * 6 + j * 2 + t] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
					}
				}
			}
		}

		for (int i = 0; i<3; i++)
		{
			for (int j = 0; j<3; j++)
			{
				for (int k = 0; k<2; k++)
				{
					for (int t = 0; t<2; t++)
					{
						if (mu_ijt[i * 6 + j * 2 + t]>0)
							mu[i * 12 + j * 4 + k * 2 + t] = mu[i * 12 + j * 4 + k * 2 + t] * n_ijt[i * 6 + j * 2 + t] / mu_ijt[i * 6 + j * 2 + t];
						else
							mu[i * 12 + j * 4 + k * 2 + t] = 0;
					}
				}
			}
		}

		//step7: for association detection
		//mu_ijkt = mu_ijkt*n_ijk/mu_ijk
		if (flag)
		{
			for (int i = 0; i<18; i++)
			{
				mu_ijk[i] = 0;
				//n_ijk[i] = 0;
			}

			for (int i = 0; i<3; i++)
			{
				for (int j = 0; j<3; j++)
				{
					for (int k = 0; k<2; k++)
					{
						for (int t = 0; t<2; t++)
						{
							mu_ijk[i * 6 + j * 2 + k] += mu[i * 12 + j * 4 + k * 2 + t];
							//n_ijk[i * 6 + j * 2 + k] += localGenoDistr[k * 18 + t * 9 + i * 3 + j];
						}
					}
				}
			}

			for (int i = 0; i<3; i++)
			{
				for (int j = 0; j<3; j++)
				{
					for (int k = 0; k<2; k++)
					{
						for (int t = 0; t<2; t++)
						{
							if (mu_ijk[i * 6 + j * 2 + k]>0)
								mu[i * 12 + j * 4 + k * 2 + t] = mu[i * 12 + j * 4 + k * 2 + t] * n_ijk[i * 6 + j * 2 + k] / mu_ijk[i * 6 + j * 2 + k];
							else
								mu[i * 12 + j * 4 + k * 2 + t] = 0;
						}
					}
				}
			}
		}

		muError = 0.0;
		for (int index = 0; index<36; index++)
		{
			muError += abs(mu[index] - mu0[index]);
		}



	}

	for (int i = 0; i<3; i++)
	{
		for (int j = 0; j<3; j++)
		{
			for (int k = 0; k<2; k++)
			{
				for (int t = 0; t<2; t++)
				{
					tmp = localGenoDistr[k * 18 + t * 9 + i * 3 + j];
					if (mu[i * 12 + j * 4 + k * 2 + t]>0)
					{
						Likelihood += tmp*log(mu[i * 12 + j * 4 + k * 2 + t]);
					}
				}
			}
		}
	}

	return Likelihood;
}

double postCorrection_2(int* localGenoDistr, bool flag)
{
	float mu[72];
	float mu0[72] = { 0 };
	float mu_ijt1[18] = { 0 };
	float mu_ijt2[18] = { 0 };
	float mu_ijk[18] = { 0 };
	float mu_t1t2[4] = { 0 };
	float mu_kt1[4] = { 0 };
	float mu_kt2[4] = { 0 };
	float mu_ij[9] = { 0 };
	float mu_ik[6] = { 0 };
	float mu_jk[6] = { 0 };
	float* vpointer[9] = { mu_ijk, mu_ijt1, mu_ijt2, mu_t1t2, mu_kt1, mu_kt2, mu_ij, mu_ik, mu_jk };


	float n_ijt1[18] = { 0 };
	float n_ijt2[18] = { 0 };
	float n_ijk[18] = { 0 };
	float n_t1t2[4] = { 0 };
	float n_kt1[4] = { 0 };
	float n_kt2[4] = { 0 };
	float n_ij[9] = { 0 };
	float n_ik[6] = { 0 };
	float n_jk[6] = { 0 };

	float* npointer[9] = { n_ijk, n_ijt1, n_ijt2, n_t1t2, n_kt1, n_kt2, n_ij, n_ik, n_jk };
	int arraySize[9] = { 18, 18, 18, 4, 4, 4, 9, 6, 6 };

	double muError = 0.0;
	double Likelihood = 0.0;
	double tmp = 0;

	int loopst, loopend;
	for (int index = 0; index < 72; index++)
	{
		mu[index] = 100;
		muError += abs(mu[index] - mu0[index]);
	}

	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						n_ijk[i * 6 + j * 2 + k] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_ijt1[i * 6 + j * 2 + t1] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_ijt2[i * 6 + j * 2 + t2] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_t1t2[t1 * 2 + t2] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_kt1[k * 2 + t1] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_kt2[k * 2 + t2] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_ij[i * 3 + j] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_ik[i * 2 + k] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						n_jk[j * 2 + k] += localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
					}
				}
			}
		}
	}


	while (muError > 1)
	{
		for (int i = 0; i < 72; i++)
		{
			mu0[i] = mu[i];
		}

		if (flag)
		{
			loopst = 0;
			loopend = 6;
		}
		else
		{
			loopst = 1;
			loopend = 9;
		}

		for (int outerloop = loopst; outerloop < loopend; outerloop++)
		{
			for (int i = 0; i < arraySize[outerloop]; i++)
			{
				vpointer[outerloop][i] = 0;
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								int index = 0;
								if (outerloop == 0)index = i * 6 + j * 2 + k;
								else if (outerloop == 1)index = i * 6 + j * 2 + t1;
								else if (outerloop == 2) index = i * 6 + j * 2 + t2;
								else if (outerloop == 3)index = t1 * 2 + t2;
								else if (outerloop == 4)index = k * 2 + t1;
								else if (outerloop == 5)index = k * 2 + t2;
								else if (outerloop == 6)index = i * 3 + j;
								else if (outerloop == 7)index = i * 2 + k;
								else if (outerloop == 8) index = j * 2 + k;

								vpointer[outerloop][index] += mu[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
							}
						}
					}
				}
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								int index = 0;
								if (outerloop == 0)index = i * 6 + j * 2 + k;
								else if (outerloop == 1)index = i * 6 + j * 2 + t1;
								else if (outerloop == 2) index = i * 6 + j * 2 + t2;
								else if (outerloop == 3)index = t1 * 2 + t2;
								else if (outerloop == 4)index = k * 2 + t1;
								else if (outerloop == 5)index = k * 2 + t2;
								else if (outerloop == 6)index = i * 3 + j;
								else if (outerloop == 7)index = i * 2 + k;
								else if (outerloop == 8) index = j * 2 + k;

								if (vpointer[outerloop][index]>0)
									mu[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j] *= (npointer[outerloop][index] / vpointer[outerloop][index]);
								else
									mu[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j] = 0;
							}
						}
					}
				}
			}

		}

		muError = 0;
		for (int i = 0; i < 72; i++)
		{
			muError += abs(mu[i] - mu0[i]);
		}

	}
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						tmp = localGenoDistr[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j];
						if (mu[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j]>0)
							Likelihood += tmp*log(mu[k * 36 + t1 * 18 + t2 * 9 + i * 3 + j]);
					}
				}
			}
		}
	}
	return Likelihood;
}

double postCorrection_3(int* localGenoDistr, bool flag)
{
	float mu[144];
	float mu0[144] = { 0 };

	float mu_ijk[18] = { 0 };	//for interacted model
	float mu_ijt1t2t3[72] = { 0 };
	float mu_kt1t2t3[16] = { 0 };
	float mu_ik[6] = { 0 };		//without interaction
	float mu_jk[6] = { 0 };		//without interaction
	float* vpointer[5] = { mu_ijk, mu_ijt1t2t3, mu_kt1t2t3, mu_ik, mu_jk };

	float n_ijk[18] = { 0 };	//for interacted model
	float n_ijt1t2t3[72] = { 0 };
	float n_kt1t2t3[16] = { 0 };
	float n_ik[6] = { 0 };		//without interaction
	float n_jk[6] = { 0 };		//without interaction
	float* npointer[5] = { n_ijk, n_ijt1t2t3, n_kt1t2t3, n_ik, n_jk };

	/*float n_ijt1[18] = { 0 };
	float n_ijt2[18] = { 0 };
	float n_ijk[18] = { 0 };
	float n_t1t2[4] = { 0 };
	float n_kt1[4] = { 0 };
	float n_kt2[4] = { 0 };
	float n_ij[9] = { 0 };
	float n_ik[6] = { 0 };
	float n_jk[6] = { 0 };

	float* npointer[9] = { n_ijk, n_ijt1, n_ijt2, n_t1t2, n_kt1, n_kt2, n_ij, n_ik, n_jk };
	*/
	int arraySize[5] = { 18, 72, 16, 6, 6 };

	double muError = 0.0;
	double Likelihood = 0.0;
	double tmp = 0;

	int loopst, loopend;
	for (int index = 0; index < 144; index++)
	{
		mu[index] = 100;
		muError += abs(mu[index] - mu0[index]);
	}

	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						for (int t3 = 0; t3 < 2; t3++)
						{
							n_ijk[i * 6 + j * 2 + k] += localGenoDistr[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j];
							n_ijt1t2t3[i * 24 + j * 8 + t1 * 4 + t2 * 2 + t3] += localGenoDistr[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j];
							n_kt1t2t3[k * 8 + t1 * 4 + t2 * 2 + t3] += localGenoDistr[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j];
							n_ik[i * 2 + k] += localGenoDistr[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j];
							n_jk[j * 2 + k] += localGenoDistr[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j];
						}
					}
				}
			}
		}
	}


	while (muError > 1)
	{
		for (int i = 0; i < 144; i++)
		{
			mu0[i] = mu[i];
		}

		if (flag)
		{
			loopst = 0;
			loopend = 3;
		}
		else
		{
			loopst = 1;
			loopend = 5;
		}

		for (int outerloop = loopst; outerloop < loopend; outerloop++)
		{
			for (int i = 0; i < arraySize[outerloop]; i++)
			{
				vpointer[outerloop][i] = 0;
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								for (int t3 = 0; t3 < 2; t3++)
								{
									int index = 0;
									if (outerloop == 0)index = i * 6 + j * 2 + k;
									else if (outerloop == 1)index = i * 24 + j * 8 + t1 * 4 + t2 * 2 + t3;
									else if (outerloop == 2) index = k * 8 + t1 * 4 + t2 * 2 + t3;
									else if (outerloop == 3)index = i * 2 + k;
									else if (outerloop == 4)index = j * 2 + k;
									vpointer[outerloop][index] += mu[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j];
								}
							}
						}
					}
				}
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								for (int t3 = 0; t3 < 2; t3++)
								{
									int index = 0;
									if (outerloop == 0)index = i * 6 + j * 2 + k;
									else if (outerloop == 1)index = i * 24 + j * 8 + t1 * 4 + t2 * 2 + t3;
									else if (outerloop == 2) index = k * 8 + t1 * 4 + t2 * 2 + t3;
									else if (outerloop == 3)index = i * 2 + k;
									else if (outerloop == 4)index = j * 2 + k;

									if (vpointer[outerloop][index]>0)
										mu[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j] *= (npointer[outerloop][index] / vpointer[outerloop][index]);
									else
										mu[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j] = 0;
								}
							}
						}
					}
				}
			}

		}

		muError = 0;
		for (int i = 0; i < 144; i++)
		{
			muError += abs(mu[i] - mu0[i]);
		}

	}
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						for (int t3 = 0; t3 < 2; t3++)
						{
							tmp = localGenoDistr[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j];
							if (mu[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j]>0)
								Likelihood += tmp*log(mu[k * 72 + t1 * 36 + t2 * 18 + t3 * 9 + i * 3 + j]);
						}
					}
				}
			}
		}
	}
	return Likelihood;
}

double postCorrection_4(int* localGenoDistr, bool flag)
{
	float mu[288];
	float mu0[288] = { 0 };

	float mu_ijk[18] = { 0 };	//for interacted model
	float mu_ijt1t2t3t4[144] = { 0 };
	float mu_kt1t2t3t4[32] = { 0 };
	float mu_ik[6] = { 0 };		//without interaction
	float mu_jk[6] = { 0 };		//without interaction
	float* vpointer[5] = { mu_ijk, mu_ijt1t2t3t4, mu_kt1t2t3t4, mu_ik, mu_jk };



	float n_ijk[18] = { 0 };	//for interacted model
	float n_ijt1t2t3t4[144] = { 0 };
	float n_kt1t2t3t4[32] = { 0 };
	float n_ik[6] = { 0 };		//without interaction
	float n_jk[6] = { 0 };		//without interaction
	float* npointer[5] = { n_ijk, n_ijt1t2t3t4, n_kt1t2t3t4, n_ik, n_jk };

	int arraySize[5] = { 18, 144, 32, 6, 6 };

	double muError = 0.0;
	double Likelihood = 0.0;
	double tmp = 0;

	int loopst, loopend;
	for (int index = 0; index < 288; index++)
	{
		mu[index] = 100;
		muError += abs(mu[index] - mu0[index]);
	}

	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						for (int t3 = 0; t3 < 2; t3++)
						{
							for (int t4 = 0; t4 < 2; t4++)
							{
								n_ijk[i * 6 + j * 2 + k] += localGenoDistr[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j];
								n_ijt1t2t3t4[i * 48 + j * 16 + t1 * 8 + t2 * 4 + t3 * 2 + t4] += localGenoDistr[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j];
								n_kt1t2t3t4[k * 16 + t1 * 8 + t2 * 4 + t3 * 2 + t4] += localGenoDistr[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j];
								n_ik[i * 2 + k] += localGenoDistr[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j];
								n_jk[j * 2 + k] += localGenoDistr[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j];
							}
						}
					}
				}
			}
		}
	}


	while (muError > 1)
	{
		for (int i = 0; i < 288; i++)
		{
			mu0[i] = mu[i];
		}

		if (flag)
		{
			loopst = 0;
			loopend = 3;
		}
		else
		{
			loopst = 1;
			loopend = 5;
		}

		for (int outerloop = loopst; outerloop < loopend; outerloop++)
		{
			for (int i = 0; i < arraySize[outerloop]; i++)
			{
				vpointer[outerloop][i] = 0;
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								for (int t3 = 0; t3 < 2; t3++)
								{
									for (int t4 = 0; t4 < 2; t4++)
									{

										int index = 0;
										if (outerloop == 0)index = i * 6 + j * 2 + k;
										else if (outerloop == 1)index = i * 48 + j * 16 + t1 * 8 + t2 * 4 + t3 * 2 + t4;
										else if (outerloop == 2) index = k * 16 + t1 * 8 + t2 * 4 + t3 * 2 + t4;
										else if (outerloop == 3)index = i * 2 + k;
										else if (outerloop == 4)index = j * 2 + k;
										vpointer[outerloop][index] += mu[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j];
									}
								}
							}
						}
					}
				}
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								for (int t3 = 0; t3 < 2; t3++)
								{
									for (int t4 = 0; t4 < 2; t4++)
									{
										int index = 0;
										if (outerloop == 0)index = i * 6 + j * 2 + k;
										else if (outerloop == 1)index = i * 48 + j * 16 + t1 * 8 + t2 * 4 + t3 * 2 + t4;
										else if (outerloop == 2) index = k * 16 + t1 * 8 + t2 * 4 + t3 * 2 + t4;
										else if (outerloop == 3)index = i * 2 + k;
										else if (outerloop == 4)index = j * 2 + k;

										if (vpointer[outerloop][index]>0)
											mu[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j] *= (npointer[outerloop][index] / vpointer[outerloop][index]);
										else
											mu[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j] = 0;
									}
								}
							}
						}
					}
				}
			}

		}

		muError = 0;
		for (int i = 0; i < 288; i++)
		{
			muError += abs(mu[i] - mu0[i]);
		}

	}
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						for (int t3 = 0; t3 < 2; t3++)
						{
							for (int t4 = 0; t4 < 2; t4++)
							{
								tmp = localGenoDistr[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j];
								if (mu[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j]>0)
									Likelihood += tmp*log(mu[k * 144 + t1 * 72 + t2 * 36 + t3 * 18 + t4 * 9 + i * 3 + j]);
							}
						}
					}
				}
			}
		}
	}
	return Likelihood;
}

 double postCorrection_5(int* localGenoDistr, bool flag)
{
	float mu[576];
	float mu0[576] = { 0 };

	float mu_ijk[18] = { 0 };	//for interacted model
	float mu_ijt1t2t3t4t5[288] = { 0 };
	float mu_kt1t2t3t4t5[64] = { 0 };
	float mu_ik[6] = { 0 };		//without interaction
	float mu_jk[6] = { 0 };		//without interaction
	float* vpointer[5] = { mu_ijk, mu_ijt1t2t3t4t5, mu_kt1t2t3t4t5, mu_ik, mu_jk };



	float n_ijk[18] = { 0 };	//for interacted model
	float n_ijt1t2t3t4t5[288] = { 0 };
	float n_kt1t2t3t4t5[64] = { 0 };
	float n_ik[6] = { 0 };		//without interaction
	float n_jk[6] = { 0 };		//without interaction
	float* npointer[5] = { n_ijk, n_ijt1t2t3t4t5, n_kt1t2t3t4t5, n_ik, n_jk };


	int arraySize[5] = { 18, 288, 64, 6, 6 };

	double muError = 0.0;
	double Likelihood = 0.0;
	double tmp = 0;

	int loopst, loopend;
	for(int index = 0; index < 576; index++)
	{
		mu[index] = 10;
		muError += abs(mu[index] - mu0[index]);
	}

	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						for (int t3 = 0; t3 < 2; t3++)
						{
							for (int t4 = 0; t4 < 2; t4++)
							{
								for (int t5 = 0; t5 < 2; t5++)
								{
									n_ijk[i * 6 + j * 2 + k] += localGenoDistr[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j];
									n_ijt1t2t3t4t5[i * 96 + j * 32 + t1 * 16 + t2 * 8 + t3 * 4 + t4 * 2 + t5] += localGenoDistr[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j];
									n_kt1t2t3t4t5[k * 32 + t1 * 16 + t2 * 8 + t3 * 4 + t4 * 2 + t5] += localGenoDistr[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j];
									n_ik[i * 2 + k] += localGenoDistr[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j];
									n_jk[j * 2 + k] += localGenoDistr[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j];
								}
							}
						}
					}
				}
			}
		}
	}


	while (muError > 1)
	{
		for (int i = 0; i < 576; i++)
		{
			mu0[i] = mu[i];
		}

		if (flag)
		{
			loopst = 0;
			loopend = 3;
		}
		else
		{
			loopst = 1;
			loopend = 5;
		}

		for (int outerloop = loopst; outerloop < loopend; outerloop++)
		{
			for (int i = 0; i < arraySize[outerloop]; i++)
			{
				vpointer[outerloop][i] = 0;
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								for (int t3 = 0; t3 < 2; t3++)
								{
									for (int t4 = 0; t4 < 2; t4++)
									{
										for (int t5 = 0; t5 < 2; t5++)
										{

											int index = 0;
											if (outerloop == 0)index = i * 6 + j * 2 + k;
											else if (outerloop == 1)index = i * 96 + j * 32 + t1 * 16 + t2 * 8 + t3 * 4 + t4 * 2 + t5;
											else if (outerloop == 2) index = k * 32 + t1 * 16 + t2 * 8 + t3 * 4 + t4 * 2 + t5;
											else if (outerloop == 3)index = i * 2 + k;
											else if (outerloop == 4)index = j * 2 + k;
											vpointer[outerloop][index] += mu[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j];
										}
									}
								}
							}
						}
					}
				}
			}

			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 2; k++)
					{
						for (int t1 = 0; t1 < 2; t1++)
						{
							for (int t2 = 0; t2 < 2; t2++)
							{
								for (int t3 = 0; t3 < 2; t3++)
								{
									for (int t4 = 0; t4 < 2; t4++)
									{
										for (int t5 = 0; t5 < 2; t5++)
										{
											int index = 0;
											if (outerloop == 0)index = i * 6 + j * 2 + k;
											else if (outerloop == 1)index = i * 96 + j * 32 + t1 * 16 + t2 * 8 + t3 * 4 + t4 * 2 + t5;
											else if (outerloop == 2) index = k * 32 + t1 * 16 + t2 * 8 + t3 * 4 + t4 * 2 + t5;
											else if (outerloop == 3)index = i * 2 + k;
											else if (outerloop == 4)index = j * 2 + k;

											if (vpointer[outerloop][index]>0)
												mu[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j] *= (npointer[outerloop][index] / vpointer[outerloop][index]);
											else
												mu[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j] = 0;
										}
									}
								}
							}
						}
					}
				}
			}

		}

		muError = 0;
		for (int i = 0; i < 576; i++)
		{
			muError += abs(mu[i] - mu0[i]);
		}

	}
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 2; k++)
			{
				for (int t1 = 0; t1 < 2; t1++)
				{
					for (int t2 = 0; t2 < 2; t2++)
					{
						for (int t3 = 0; t3 < 2; t3++)
						{
							for (int t4 = 0; t4 < 2; t4++)
							{
								for (int t5 = 0; t5 < 2; t5++)
								{
									tmp = localGenoDistr[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j];
									if (mu[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j]>0)
										Likelihood += tmp*log(mu[k * 288 + t1 * 144 + t2 * 72 + t3 * 36 + t4 * 18 + t5 * 9 + i * 3 + j]);
								}
							}
						}
					}
				}
			}
		}
	}
	return Likelihood;
}
