using System;

namespace Science.Biology.PopulationGenetics
{
	/// <summary>
	/// PopulationEquation
	/// </summary>
	public class PopulationEquation
	{
		private double[] adultXX;
		private double[] adultXY;
		private double[] adultPhenoXX;
		private double[] adultPhenoXY;
		private double[] birthXX;
		private double[] birthXY;
		private double[,] parents;
		private double[,] effectiveparents;
		private double[,] breed;
		private double[,] Gamma;
		private double[,] GammaPheno;
		private int[] mapFromGenoToPhenoXX;
		private int[] mapFromGenoToPhenoXY;
		private double[,] mutationXX;
		private double[,] mutationXY;
		private double[,,] reproductioncoefficientXX;
		private double[,,] reproductioncoefficientXY;
		private double[] disadvantagefactorXX;
		private double[] disadvantagefactorXY;
		private bool matingtype;
		private bool phenoMating;
		private bool hierarchicalMating;

		public double[] AdultFemalePopulation
		{
			get{return adultXX;}
		}
		public double[] AdultMalePopulation
		{
			get{return adultXY;}
		}
		public double[] AdultPhenoFemalePopulation
		{
			get{return adultPhenoXX;}
		}
		public double[] AdultPhenoMalePopulation
		{
			get{return adultPhenoXY;}
		}
		public double[] BirthFemalePopulation
		{
			get{return birthXX;}
		}
		public double[] BirthMalePopulation
		{
			get{return birthXY;}
		}
		public double[,] ParentPopulation
		{
			get{return parents;}
		}
		public double[,] MatingFactor
		{
			get{return Gamma;}
		}

		public PopulationEquation(Mating m, Mutation muta, Reproduction r, Selection s)
		{
			adultXX = new double[muta.NumberOfAlleleFemale];
			adultXY = new double[muta.NumberOfAlleleMale];
			birthXX = new double[muta.NumberOfAlleleFemale];
			birthXY = new double[muta.NumberOfAlleleMale];
			parents = new double[muta.NumberOfAlleleFemale, muta.NumberOfAlleleMale];
            effectiveparents = new double[muta.NumberOfAlleleFemale, muta.NumberOfAlleleMale];

			Gamma = new double[muta.NumberOfAlleleFemale, muta.NumberOfAlleleMale];
			for(int i = 0; i < muta.NumberOfAlleleFemale ; i++) 
			{
				for(int j = 0; j < muta.NumberOfAlleleMale ; j++) 
				{
					Gamma[i,j] = 1.0;
				}	
			}
            phenoMating = m.PhenotypeQ;
			matingtype = m.ObjectiveQ;
			hierarchicalMating = m.HierarchicalQ;
            if (matingtype) breed = m.BreedParameter;
            if (phenoMating)
            {
                adultPhenoXX = new double[m.NumberOfPhenoAlleleFemale];
                adultPhenoXY = new double[m.NumberOfPhenoAlleleMale];
                GammaPheno = new double[m.NumberOfPhenoAlleleFemale, m.NumberOfPhenoAlleleMale];
                mapFromGenoToPhenoXX = m.MapFromGenoToPhenoFemale;
                mapFromGenoToPhenoXY = m.MapFromGenoToPhenoMale;		

            }

			mutationXX = muta.MutationFemale;
			mutationXY = muta.MutationMale;	

			reproductioncoefficientXX = r.CoefficientXX;
			reproductioncoefficientXY = r.CoefficientXY;

			disadvantagefactorXX = s.DisadvantageXX;
			disadvantagefactorXY = s.DisadvantageXY;

			InitialRun();
     		Normalization();		
		}
	
		private void InitialRun()
		{
			if(hierarchicalMating)
			{
				Random ran = new Random();
				for(int i = 0; i < adultXX.Length; i++) 
				{
					adultXX[i] = (double) ran.Next(0,100);
					adultXY[i] = adultXX[i]+0.01;
				}
			}
			else
			{
				Random ran = new Random();
				for(int i = 0; i < adultXX.Length; i++) adultXX[i] = (double) ran.Next(0,100);
				for(int i = 0; i < adultXY.Length; i++) adultXY[i] = (double) ran.Next(0,100);
			}
		}

		private void Normalization()
		{
			double sum = 0.0;
			for(int i = 0; i < adultXX.Length; i++) sum += adultXX[i];
			for(int i = 0; i < adultXX.Length; i++) adultXX[i] /= sum;
			sum = 0.0;
			for(int i = 0; i < adultXY.Length; i++) sum += adultXY[i];
			for(int i = 0; i < adultXY.Length; i++) adultXY[i] /= sum;
		}

		private double[] initialPopulationXX, initialPopulationXY;
		public double[] InitialAdultMalePopulation
		{
			set{initialPopulationXY = value;}
		}
		public double[] InitialAdultFemalePopulation
		{
			set{initialPopulationXX = value;}
		}
		public void SetInitialPopulation()
		{
			for(int i = 0; i < adultXX.Length; i++) adultXX[i] = initialPopulationXX[i];
			for(int i = 0; i < adultXY.Length; i++) adultXY[i] = initialPopulationXY[i];
		}

		private int howManyTimes;
		public int HowManyTimes
		{
			set{howManyTimes = value;}
		}
		public void Compute()  // after set
		{
            for (int i = 0; i < howManyTimes; i++)
            {
                Normalization();
                MatingStep();
                MutationStep();
                ReproductionStep();
                SelectionStep();
            }          		
		}  

		private void MatingStep()
		{
			if(!matingtype & !hierarchicalMating)   // random
			{
				for(int i = 0; i < adultXX.Length; i++)
					for(int j = 0; j < adultXY.Length; j++) parents[i,j] = adultXX[i]*adultXY[j];

			}
			else   // non-random
			{
				if(!phenoMating)  // genotype
				{
					if(!hierarchicalMating)  // breedfunctiontype
					{						
						LinearOptimization(breed,adultXX,adultXY,Gamma);       // Find Gamma[,] 
						for(int i = 0; i < adultXX.Length; i++)
							for(int j = 0; j < adultXY.Length; j++) parents[i,j] = adultXX[i]*Gamma[i,j]*adultXY[j];
					}
					else  // hierarchicaltype
					{
						for(int i = 0; i < adultXX.Length; i++)
							for(int j = 0; j < adultXY.Length; j++) parents[i,j] = 0.0;
						parents[0,0] = adultXX[0];
						for(int i = 0; i < adultXX.Length - 1; i++) 
						{
							parents[i+1,i] = adultXY[i] - parents[i,i];
							parents[i+1,i+1] = adultXX[i+1] - parents[i+1,i];
						}
					}
				}
				else   // phenotype
				{
					for(int i = 0; i < adultPhenoXX.Length; i++)
					{						
						adultPhenoXX[i] = 0.0;
					}
					for(int i = 0; i < adultXX.Length; i++)
					{						
						adultPhenoXX[mapFromGenoToPhenoXX[i]] += adultXX[i];
					}
					for(int i = 0; i < adultPhenoXY.Length; i++)
					{						
						adultPhenoXY[i] = 0.0;
					}
					for(int i = 0; i < adultXY.Length; i++)
					{						
						adultPhenoXY[mapFromGenoToPhenoXY[i]] += adultXY[i];
					}
					if(!hierarchicalMating)  // breedfunctiontype
					{						
						LinearOptimization(breed,adultPhenoXX,adultPhenoXY,GammaPheno);
						for(int i = 0; i < adultXX.Length; i++)
							for(int j = 0; j < adultXY.Length; j++) 
							{
								Gamma[i,j]=GammaPheno[mapFromGenoToPhenoXX[i],mapFromGenoToPhenoXY[j]];
								parents[i,j] = adultXX[i]*Gamma[i,j]*adultXY[j];
							}
					}
					else  // hierarchicaltype
					{
						for(int i = 0; i < adultPhenoXX.Length; i++)
							for(int j = 0; j < adultPhenoXY.Length; j++) GammaPheno[i,j] = 0.0;
						GammaPheno[0,0] = 1.0/adultPhenoXY[0];
						for(int i = 0; i < adultPhenoXX.Length - 1; i++) 
						{
							GammaPheno[i+1,i] = (1.0 - adultPhenoXX[i]*GammaPheno[i,i])/adultPhenoXX[i+1];
							GammaPheno[i+1,i+1] = (1.0 - GammaPheno[i+1,i]*adultPhenoXY[i])/adultPhenoXY[i+1];
						}
						for(int i = 0; i < adultXX.Length; i++)
							for(int j = 0; j < adultXY.Length; j++) 
							{
								Gamma[i,j]=GammaPheno[mapFromGenoToPhenoXX[i],mapFromGenoToPhenoXY[j]];
								parents[i,j] = adultXX[i]*Gamma[i,j]*adultXY[j];
							}
					}
				}
			}	
		}
		
		private void LinearOptimization(double[,] bp, double[] pf, double[] pm, double[,] mf)
		{
			double[,] a;
			int[] izrov, iposv;
			
			if(adultXX.Length == adultXY.Length)  // auto chromosome
			{
				a = new Double[2+pf.Length+pf.Length*(pf.Length-1)/2,1+pf.Length*pm.Length];
				izrov = new Int32[pf.Length*pm.Length];
				iposv = new Int32[pf.Length+pf.Length*(pf.Length-1)/2];

				int kk = 0;
				for(int i = 0; i < pf.Length ; i++) 
				{
					for(int j = 0; j < pm.Length ; j++) 
					{
						a[0,1+kk] = bp[i,j]*pf[i]*pm[j];
						kk++;
					}	
				}

     			for(int i = 0; i < pf.Length; i++)
				{
					a[1+i,0] = 1.0;
					for(int j = 0; j < pm.Length ; j++) a[1+i,1+pm.Length*i+j] = -pm[j];			
				}

				int k = 0;
				for(int i = 0; i < pf.Length ; i++) 
				{
					for(int j = 0; j < i ; j++) 
					{
						a[1+pf.Length+k,1+pm.Length*i+j] = -1.0;
						a[1+pf.Length+k,1+pm.Length*j+i] = 1.0;
						k++;
					}	
				}
      
				simplx(a,pf.Length+pf.Length*(pf.Length-1)/2,pf.Length*pm.Length,0,0,pf.Length+pf.Length*(pf.Length-1)/2,izrov,iposv);
		
				for(int i = 0; i < pf.Length ; i++) 
				{
					for(int j = 0; j < pm.Length ; j++) 
					{
						mf[i,j] = 0.0;
					}	
				}

				for(int i = 0; i < iposv.Length; i++) 
				{
					int v = iposv[i];
					int jj, ii;
					if(v <= pf.Length*pm.Length)
					{
						jj = (v-1) % pm.Length;
						ii = (v-1-jj)/pm.Length;
						mf[ii,jj] = a[1+i,0];
					}
				}		
			}
			else         // sex-linked
			{
				a = new Double[2+pf.Length+pm.Length,1+pf.Length*pm.Length];
				izrov = new int[pf.Length*pm.Length];
				iposv = new int[pf.Length+pm.Length];
				int kk = 0;
				for(int i = 0; i < pf.Length ; i++) 
				{
					for(int j = 0; j < pm.Length ; j++) 
					{
						a[0,1+kk] = bp[i,j]*pf[i]*pm[j];
						kk++;
					}	
				}
				for(int i = 0; i < pf.Length; i++)
				{
					a[1+i,0] = 1.0;
					for(int j = 0; j < pm.Length; j++)
						a[1+i,1+pm.Length*i+j] = -pm[j];
				}
				for(int i = 0; i < pm.Length; i++)
				{
					a[1+pf.Length+i,0] = 1.0;
					for(int j = 0; j < pf.Length; j++)
						a[1+pf.Length+i,1+i+pm.Length*j] = -pf[j];
				}

				simplx(a,pf.Length+pm.Length,pf.Length*pm.Length,0,0,pf.Length+pm.Length,izrov,iposv);

				for(int i = 0; i < pf.Length ; i++) 
				{
					for(int j = 0; j < pm.Length ; j++) 
					{
						mf[i,j] = 0.0;
					}	
				}

				for(int i = 0; i < iposv.Length; i++) 
				{
					int k = iposv[i];
					if(k <= pf.Length*pm.Length)
					{
						int jj = (k-1) % pm.Length;
						int ii = ((k-1) - ((k-1) % pm.Length))/pm.Length;
						mf[ii,jj] = a[1+i,0];
					}
				}
			}
		}	
	
		private void MutationStep()
		{
			double sum;
			for(int i = 0; i < adultXX.Length; i++)
				for(int j = 0; j < adultXY.Length; j++)
				{
					sum = 0.0;
					for(int k = 0; k < adultXX.Length; k++)
						for(int l = 0; l < adultXY.Length; l++)
							sum += parents[k,l]*mutationXX[k,i]*mutationXY[l,j];
					effectiveparents[i,j] = sum;
				}
		}

		private void ReproductionStep()
		{
			double sum;
			for(int i = 0; i < adultXX.Length; i++)
			{  
				sum = 0.0;
				for(int j = 0; j < adultXX.Length; j++)
					for(int k = 0; k < adultXY.Length; k++)
						sum  += effectiveparents[j,k]*reproductioncoefficientXX[j,k,i];
				birthXX[i] = sum;
			}	
			for(int i = 0; i < adultXY.Length; i++)
			{  
				sum = 0.0;
				for(int j = 0; j < adultXX.Length; j++)
					for(int k = 0; k < adultXY.Length; k++)
						sum  += effectiveparents[j,k]*reproductioncoefficientXY[j,k,i];
				birthXY[i] = sum;
			}	
		}

		private void SelectionStep()
		{
			double sum = 0.0;
			for(int k = 0; k < adultXX.Length; k++)
				sum  += birthXX[k]*disadvantagefactorXX[k];
			for(int i = 0; i < adultXX.Length; i++)
			    adultXX[i] = birthXX[i]*(1.0 - disadvantagefactorXX[i])/(1.0 - sum);
			sum = 0.0;
			for(int k = 0; k < adultXY.Length; k++)
				sum  += birthXY[k]*disadvantagefactorXY[k];
			for(int i = 0; i < adultXY.Length; i++)
				adultXY[i] = birthXY[i]*(1.0 - disadvantagefactorXY[i])/(1.0 - sum);
		}

		public string InfantMortality()
		{
			string result;
            double sum = 0.0;
			for(int k = 0; k < adultXX.Length; k++)
				if(disadvantagefactorXX[k] > 0.5)  sum  += birthXX[k];
			result = "Female\r\n";
			result += "Normal vs Defect = "+Convert.ToString(1.0 - sum)+" : "+Convert.ToString(sum)+"\r\n";
			sum = 0.0;
			for(int k = 0; k < adultXY.Length; k++)
				if(disadvantagefactorXY[k] > 0.5)  sum  += birthXY[k];
			result += "Male\r\n";
			result += "Normal vs Defect = "+Convert.ToString(1.0 - sum)+" : "+Convert.ToString(sum)+"\r\n";
            return result;
		}

// numerical recipes source code
		private double EPS = 1.0e-6;

		private int icase;
		public int SimplexMethodResult
		{
			get{return icase;}
		}
		private int kp = 0;
		private double bmax = 0.0;
		private int ip = 0;

		private void simplx(double[,] a, int m, int n, int m1, int m2, int m3, 
			int[] izrov, int[] iposv)
		{
			int i,iis,k,kh,nl1;
			int[] l1 = new Int32[n+1];
			int[] l3 = new Int32[m];
			double q1;

            if (m != (m1 + m2 + m3)) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad input constraint counts in simplx");
                } 
			nl1=n;
			for (k=1;k<=n;k++) l1[k-1]=izrov[k-1]=k;
			for (i=1;i<=m;i++) 
			{
                if (a[i, 0] < 0.0) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Bad input tableau in simplx");
                    }
				iposv[i-1]=n+i;
			}
			if ((m2+m3) != 0) 
			{
				for (i=1;i<=m2;i++) l3[i-1]=1;
				for (k=1;k<=(n+1);k++) 
				{
					q1=0.0;
					for (i=m1+1;i<=m;i++) q1 += a[i,k-1];
					a[m+1,k-1] = -q1;
				}
				for (;;) 
				{
					simp1(a,m+1,l1,nl1,0);
					if (bmax <= EPS && a[m+1,0] < -EPS) 
					{
						icase = -1;
						return;
					} 
					else if (bmax <= EPS && a[m+1,0] <= EPS) 
					{
						for (ip=m1+m2+1;ip<=m;ip++) 
						{
							if (iposv[ip-1] == (ip+n)) 
							{
								simp1(a,ip,l1,nl1,1);
								if (bmax > EPS)
								goto one;
							}
						}
						for (i=m1+1;i<=m1+m2;i++)
							if (l3[i-m1-1] == 1)
								for (k=1;k<=n+1;k++)
									a[i,k-1] = -a[i,k-1];
						break;
					}
				    simp2(a,m,n,kp);
					if (ip == 0) 
					{
						icase = -1;
						return;
					}
				one:	simp3(a,m+1,n,ip,kp);
					if (iposv[ip-1] >= (n+m1+m2+1)) 
					{
						for (k=1;k<=nl1;k++)
							if (l1[k-1] == kp) break;
						--nl1;
						for (iis=k;iis<=nl1;iis++) l1[iis-1]=l1[iis+1-1];
					} 
					else 
					{
						kh=iposv[ip-1]-m1-n;
						if ((kh >= 1) && (l3[kh-1] != 0)) 
						{
							l3[kh-1]=0;
							++a[m+1,kp];
							for (i=1;i<=m+2;i++)
								a[i-1,kp] = -a[i-1,kp];
						}
					}
					iis=izrov[kp-1];
					izrov[kp-1]=iposv[ip-1];
					iposv[ip-1]=iis;
				}
			}
			for (;;) 
			{
				simp1(a,0,l1,nl1,0);
				if (bmax <= EPS) 
				{
					icase=0;
					return;
				}
				simp2(a,m,n,kp);
				if (ip == 0) 
				{
					icase=1;
					return;
				}
				simp3(a,m,n,ip,kp);
				iis=izrov[kp-1];
				izrov[kp-1]=iposv[ip-1];
				iposv[ip-1]=iis;
			}
		}  

	
		private void simp1(double[,] a, int mm, int[] ll, int nll, int iabf)
		{
			int k;
			double test;

			if (nll <= 0)
				bmax=0.0;
			else 
			{
				kp=ll[0];
				bmax=a[mm,kp];
				for (k=1;k<nll;k++) 
				{
					if (iabf == 0)
						test=a[mm,ll[k]]-bmax;
					else
						test=Math.Abs(a[mm,ll[k]])-Math.Abs(bmax);
					if (test > 0.0) 
					{
						bmax=a[mm,ll[k]];
						kp=ll[k];
					}
				}
			}
		} 


	
		private void simp2(double[,] a, int m, int n, int kp)
		{
			int k,i;
			double qp=0.0,q0=0.0,q,q1;

			ip=0;
			for (i=1;i<=m;i++)
				if (a[i,kp] < -EPS) break;
			if (i>m) return;
			q1 = -a[i,0]/a[i,kp];
			ip=i;
			for (i=ip+1;i<=m;i++) 
			{
				if (a[i,kp] < -EPS) 
				{
					q = -a[i,0]/a[i,kp];
					if (q < q1) 
					{
						ip=i;
						q1=q;
					} 
					else if (q == q1) 
					{
						for (k=1;k<=n;k++) 
						{
							qp = -a[ip,k]/a[ip,kp];
							q0 = -a[i,k]/a[i,kp];
							if (q0 != qp) break;
						}
						if (q0 < qp) ip=i;
					}
				}
			}
		}   


		private void simp3(double[,] a, int i1, int k1, int ip, int kp)
		{
			int kk,ii;
			double piv;

			piv=1.0/a[ip,kp];
			for (ii=1;ii<=i1+1;ii++)
				if (ii-1 != ip) 
				{
					a[ii-1,kp] *= piv;
					for (kk=1;kk<=k1+1;kk++)
						if (kk-1 != kp)
							a[ii-1,kk-1] -= a[ip,kk-1]*a[ii-1,kp];
				}
			for (kk=1;kk<=k1+1;kk++)
				if (kk-1 != kp) a[ip,kk-1] *= -piv;
			a[ip,kp]=piv;
		}  

	}
}