using System;

namespace Science.Biology.PopulationGenetics
{
	/// <summary>
	/// Reproduction
	/// </summary>
	public class Reproduction
	{
		private double[,,] coefficientXX;
		private double[,,] coefficientXY;
		public double[,,] CoefficientXX
		{
			get{return coefficientXX;}
		}
		public double[,,] CoefficientXY
		{
			get{return coefficientXY;}
		}
		private GenotypeFemale ff;
		private GenotypeMale mm;
		private SocialGene[] sg;
		public Reproduction(GenotypeFemale fg, GenotypeMale mg) //, SocialGene[] socialgene)
		{
			ff = fg;
			mm = mg;
			sg = fg.SocialGenes;
			coefficientXX = new double[fg.NumberOfAllele,mg.NumberOfAllele,fg.NumberOfAllele];
			coefficientXY = new double[fg.NumberOfAllele,mg.NumberOfAllele,mg.NumberOfAllele];
			
			for(int i = 0; i < fg.NumberOfAllele; i++) 
			{
				for(int j = 0; j < mg.NumberOfAllele; j++)
				{
					for(int k = 0; k < fg.NumberOfAllele; k++)
						coefficientXX[i,j,k] = FindXX(i,j,k);
					for(int k = 0; k < mg.NumberOfAllele; k++)
						coefficientXY[i,j,k] = FindXY(i,j,k);
				}
			}
		}

		private double FindXX(int a, int b, int c)
		{
			int v,w,p,q,s,t;
			double sum;
			double prod = 1.0;
			for(int k = 0; k < mm.ConfigurationDetail.GetLength(1); k++)
			{
				v = ff.ConfigurationDetail[a,k,0];
				w = ff.ConfigurationDetail[a,k,1];
				p = mm.ConfigurationDetail[b,k,0];
				q = mm.ConfigurationDetail[b,k,1];
				s = ff.ConfigurationDetail[c,k,0];
				t = ff.ConfigurationDetail[c,k,1];
				sum = 0.0;
				if( s == Math.Min(v,p) && t == Math.Max(v,p)) sum += 0.25;
				if( s == Math.Min(v,q) && t == Math.Max(v,q)) sum += 0.25;
				if( s == Math.Min(w,p) && t == Math.Max(w,p)) sum += 0.25;
				if( s == Math.Min(w,q) && t == Math.Max(w,q)) sum += 0.25;
				prod *= sum;
			}
			v = ff.ConfigurationDetail[a,ff.ConfigurationDetail.GetLength(1)-1,0];
			w = ff.ConfigurationDetail[a,ff.ConfigurationDetail.GetLength(1)-1,1];
			p = mm.Configuration[b,mm.ConfigurationDetail.GetLength(1)];
			s = ff.ConfigurationDetail[c,ff.ConfigurationDetail.GetLength(1)-1,0];
			t = ff.ConfigurationDetail[c,ff.ConfigurationDetail.GetLength(1)-1,1];
			sum = 0.0;
			if( s == Math.Min(v,p) && t == Math.Max(v,p)) sum += 0.5;
			if( s == Math.Min(w,p) && t == Math.Max(w,p)) sum += 0.5;
			prod *= sum;
			sum = 0.0;
			v = ff.Configuration[a,ff.ConfigurationDetail.GetLength(1)];
			s = ff.Configuration[c,ff.ConfigurationDetail.GetLength(1)];
	        if(s == v) sum = 1.0;
		    prod *= sum;
			for(int ii = 0; ii < sg.Length; ii++)
			{
				v = ff.Configuration[a,ff.ConfigurationDetail.GetLength(1)+1+ii];
				p = mm.Configuration[b,mm.ConfigurationDetail.GetLength(1)+3+ii];
				s = ff.Configuration[c,ff.ConfigurationDetail.GetLength(1)+1+ii];
				prod *= sg[ii].ReproductionXX[v,p,s];
			}
			return prod;
        }	

		private double FindXY(int a, int b, int c)
		{
			int v,w,p,q,s,t;
			double sum;
			double prod = 1.0;
			for(int k = 0; k < mm.ConfigurationDetail.GetLength(1); k++)
			{
				v = ff.ConfigurationDetail[a,k,0];
				w = ff.ConfigurationDetail[a,k,1];
				p = mm.ConfigurationDetail[b,k,0];
				q = mm.ConfigurationDetail[b,k,1];
				s = mm.ConfigurationDetail[c,k,0];
				t = mm.ConfigurationDetail[c,k,1];
				sum =0.0;
				if( s == Math.Min(v,p) && t == Math.Max(v,p)) sum += 0.25;
				if( s == Math.Min(v,q) && t == Math.Max(v,q)) sum += 0.25;
				if( s == Math.Min(w,p) && t == Math.Max(w,p)) sum += 0.25;
				if( s == Math.Min(w,q) && t == Math.Max(w,q)) sum += 0.25;
				prod *= sum;
			}
			v = ff.ConfigurationDetail[a,ff.ConfigurationDetail.GetLength(1)-1,0];
			w = ff.ConfigurationDetail[a,ff.ConfigurationDetail.GetLength(1)-1,1];
			q = mm.Configuration[b,mm.ConfigurationDetail.GetLength(1)+1];
			s = mm.Configuration[c,mm.ConfigurationDetail.GetLength(1)];
			t = mm.Configuration[c,mm.ConfigurationDetail.GetLength(1)+1];
			sum =0.0;
			if( s == v && t == q) sum += 0.5;
			if( s == w && t == q) sum += 0.5;
			prod *= sum;
			sum = 0.0;
			v = ff.Configuration[a,ff.ConfigurationDetail.GetLength(1)];
			s = mm.Configuration[c,mm.ConfigurationDetail.GetLength(1)+2];
			if(s == v) sum = 1.0;
			prod *= sum;
			for(int ii = 0; ii < sg.Length; ii++)
			{
				v = ff.Configuration[a,ff.ConfigurationDetail.GetLength(1)+1+ii];
				p = mm.Configuration[b,mm.ConfigurationDetail.GetLength(1)+3+ii];
				s = mm.Configuration[c,mm.ConfigurationDetail.GetLength(1)+3+ii];
				prod *= sg[ii].ReproductionXY[v,p,s];
			}
			return prod;
		}	
	
		private int whichFemale, whichMale, whichFemaleBaby, whichMaleBaby;
		private double modifiedvalue;
		public int WhichFemale
		{
			set{whichFemale = value;}
		}
		public int WhichMale
		{
			set{whichMale = value;}
		}
		public int WhichFemaleBaby
		{
			set{whichFemaleBaby = value;}
		}
		public int WhichMaleBaby
		{
			set{whichMaleBaby = value;}
		}
		public double ModifiedValue
		{
			set{modifiedvalue = value;}
		}
		public void SetCoefficientXX()  // after set
		{
			coefficientXX[whichFemale,whichMale,whichFemaleBaby] = modifiedvalue;
		}
		public void SetCoefficientXY()  // after set
		{
			coefficientXY[whichFemale,whichMale,whichMaleBaby] = modifiedvalue;
		}
	}
}