////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Performance evaluation of an improved harmony search algorithm for numerical optimization : Melody Search (MS)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//S.M. Ashrafi , A.B.Dariane////////////////////////////////////////////////////////////////////////////////////
// Programmed by : Ošuz Emrah Turgut 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

abstract class f_xj
{
abstract double func(double x[]);
}

public class melody_search
{
int D;// Dimension(Number of pitch)	
int PMN;//number of player memory
int PMS;//player memory size
int Maxiter1;//maximum of iteration for  first phase
int Maxiter2;//maximum number of iteration for second phase
double PMCR;//player memory considering rate
double PARmax;// Maximum value of pitch adjustment rate
double PARmin;// Minimum value of pitch adjustment rate
double PAR;//Pitch adjustment rate(chenges within every iteration) 
double Xnew[][];//Holds new pitches for stage 1
double Xnew2[];//Holds new pitches for stage 2
double[][] fitness;// Holds fitness values for stage 1
double[][] fitness2;// Holds fitness value for stage 2
double[] fitness3;
double[][][] MM;//Melody memory first stage
double[][] MM2;//Melody memory second stage
double[] Lower1;// Lower bound of first stage 
double[] Upper1;// Upper bound of first stage
double[] Lower2;// Lower bound of second stage (determined by first stage)
double[] Upper2;// Upper bound of second stage (determined by first stage)
double[] maxval;
int[] maxindex;
double[] minval;
int[] minindex; 
double[] maxval2;
int[] maxindex2;
double[] minval2;
int[] minindex2; 
double[] bw;//bandwith
double[][] worst;
double[][] best;
double[][] best2;
double[] best3;
double[] worst3;
int L;//random integer between 1 and PMS
int H;//random integer between 1 and D
double min2;
int mindex2;
double max2;
int maxdex2;
f_xj ff;

        public melody_search(f_xj iff,int iPMN,int iPMS,int iMaxiter1,int iMaxiter2,double iPMCR,double[] iLower,double[] iUpper)
        {
	    ff=iff;   
	    PMN=iPMN;    
	    PMS=iPMS;    
	    Maxiter1=iMaxiter1;    
	    Maxiter2=iMaxiter2;
	    PMCR=iPMCR;
	    Lower1=iLower;    
	    Upper1=iUpper;
	    D=Upper1.length;
	    Lower2=new double[D];
	    Upper2=new double[D];
	    maxval=new double[PMN];
	    maxindex=new int[PMN];
	    minval=new double[PMN];
	    minindex=new int[PMN];
	    maxval2=new double[PMN];
	    maxindex2=new int[PMN];
	    minval2=new double[PMN];
	    minindex2=new int[PMN];
	    fitness=new double[PMN][PMS];
	    fitness2=new double[PMN][PMS];
	    fitness3=new double[PMN];
	    worst=new double[PMN][D];
	    best=new double[PMN][D];
	    best2=new double[PMN][D];
	    best3=new double[D];
	    worst3=new double[D];
	    MM=new double[PMN][PMS][D];
	    MM2=new double[PMN][D];
	    Xnew=new double[PMN][D];
	    Xnew2=new double[D];
	    PARmax=2.0;
	    PARmin=0.0;
	    
	    }

	    void initialize()
	    {
		// Improvisation a new melody    
		for(int i=0;i<PMN;i++)
		{for(int j=0;j<PMS;j++)
		 {for(int k=0;k<D;k++)
		   {MM[i][j][k]=Lower1[k]+(Math.random()*(Upper1[k]-Lower1[k]));}}}    
		   
		   //Calculate fitness of each player
		   for(int i=0;i<PMN;i++)
		   {for(int j=0;j<PMS;j++)
		    {fitness[i][j]=ff.func(MM[i][j]);}}
		    
		   // Calculate best and worst melody for each player
		     for(int i=0;i<PMN;i++)
		     {
			  double dd[]=getmaxval_index(fitness[i]);
			  maxval[i]=dd[0];
			  maxindex[i]=(int)dd[1];      
			 }
			 
			 for(int i=0;i<PMN;i++)
		     {
			  double dd2[]=getminval_index(fitness[i]);
			  minval[i]=dd2[0];
			  minindex[i]=(int)dd2[1];      
			 }
			 for(int i=0;i<PMN;i++)
			 {
			    for(int k=0;k<D;k++)
			    {
				best[i][k]=MM[i][minindex[i]][k];    
				}
			 }
			 for(int i=0;i<PMN;i++)
			 {
			    for(int k=0;k<D;k++)
			    {
				worst[i][k]=MM[i][maxindex[i]][k];    
				}
			 }
		}
	    
		
		double[][] solution()
		{
		initialize();
		
		/// FIRST STAGE /// Determine best melody for each music player
			
		int iter1=1;	
			while(iter1<Maxiter1)
			{
			     for(int i=0;i<PMN;i++)
			     {
				     for(int k=0;k<D;k++)
				     {
					    if(Math.random()<=PMCR) // memory consideration   
			            {
	                        L=(int)(((double)PMS)*Math.random());			             
				            H=(int)(((double)D)*Math.random());
	                        PAR=PARmax-((PARmax-PARmin)*((double)iter1)/((double)Maxiter1));
				            if((iter1%2)==1)
				            {
					        Xnew[i][k]=MM[i][L][k]+(((Math.random()-0.5)*2.0)*((Upper1[k]-Lower1[k])/200.0));   
					        }
				            else
				            {
					        Xnew[i][k]=MM[i][L][H]+(((Math.random()-0.5)*2.0)*((Upper1[k]-Lower1[k])/200.0));    
					        }
				            if(Math.random()<=PAR) // pitch adjustment
				            {
					        Xnew[i][k]=best[i][k];  
					        }
				            
				        }		     
					    else // randomization
					    {
						Xnew[i][k]=Lower1[k]+(Math.random()*(Upper1[k]-Lower1[k]));    
						} 
		             }     
				 }	
				
				 for(int i=0;i<PMN;i++)   // Xnew is better than worst , exclude worst 
				 {if(ff.func(Xnew[i])<ff.func(worst[i]))
				     {for(int k=0;k<D;k++)
					    {worst[i][k]=Xnew[i][k];}
					 }
				 }
				 
				 // Update  worst and best melody 
				 for(int i=0;i<PMN;i++)
		            {for(int j=0;j<PMS;j++)
		               {fitness[i][j]=ff.func(MM[i][j]);}}
		    
		   
		         for(int i=0;i<PMN;i++)
		         {
			     double dd3[]=getmaxval_index(fitness[i]);
			     maxval[i]=dd3[0];
			     maxindex[i]=(int)dd3[1];      
			     }
			 
			     for(int i=0;i<PMN;i++)
		         {
			     double dd4[]=getminval_index(fitness[i]);
			     minval[i]=dd4[0];
			     minindex[i]=(int)dd4[1];      
			     }
				 
			     for(int i=0;i<PMN;i++)
			     {
			        for(int k=0;k<D;k++)
			        {
				     best[i][k]=MM[i][minindex[i]][k];    
				    }
			     }
				
				 for(int i=0;i<PMN;i++)
			     {
			        for(int k=0;k<D;k++)
			        {
		             MM[i][maxindex[i]][k]=worst[i][k]; // worst is now  new maxindexed MM 
				    }
			     }
				 
			     iter1++;
			}
			
			
			/// SECOND STAGE /// Determine best harmony from selected melodies
			
			
			
			// Determine new upper and lower bound 
			for(int i=0;i<PMN;i++)
			{
			   for(int j=0;j<PMS;j++)	
			    {fitness2[i][j]=ff.func(MM[i][j]);}	
		    }
			
			
			 
			 for(int i=0;i<PMN;i++)
		     {
			  double dd6[]=getminval_index(fitness2[i]);
			  minval2[i]=dd6[0];
			  minindex2[i]=(int)dd6[1];      
			 }
			 for(int i=0;i<PMN;i++)
			 {
			    for(int k=0;k<D;k++)
			    {
				best2[i][k]=MM[i][minindex[i]][k];    
				}
			 }
			 
			 for(int i=0;i<PMN;i++)
			 {
			     for(int k=0;k<D;k++)
			     {
				     MM2[i][k]=best[i][k];
				 }	 
			 }
			 
			 // MM2 now is new harmony matrix
			 
			 
			 // Procedure for determining lower and upper bounds for every column
			 double min=0.0;
			 for(int k=0;k<D;k++)
			 {
				 min=MM2[0][k];
			     for(int i=0;i<PMN;i++)
			     {
				  if(MM2[i][k]<min)
				  {min=MM2[i][k];}    
				 }
				 Lower2[k]=min;	 
			 }
			 
			 
			 
			 double max=0.0;
			 for(int k=0;k<D;k++)
			 {
				 max=MM2[0][k];
			     for(int i=0;i<PMN;i++)
			     {
				  if(MM2[i][k]>max)
				  {max=MM2[i][k];}    
				 }
				 Upper2[k]=max;	 
			 }
			
			
			// determine best and worst vector in MM2
			for(int i=0;i<PMS;i++)
			{fitness3[i]=ff.func(MM2[i]);}
			
			double dd7[]=getminval_index(fitness3);
			min2=dd7[0];
			mindex2=(int)dd7[1];
			
			double dd8[]=getmaxval_index(fitness3);
			max2=dd8[0];
			maxdex2=(int)dd8[1];
			
			for(int i=0;i<D;i++)
			{best3[i]=MM2[mindex2][i];
			 worst3[i]=MM2[maxdex2][i];}
			
			 
			         
			
			int iter2=1;
			while(iter2<Maxiter2)
			{
				
				   for(int k=0;k<D;k++)
				   {
					   if(Math.random()<PMCR)//Memory consideration
					   {
						 L=(int)(((double)PMS)*Math.random());			             
				         H=(int)(((double)D)*Math.random());  
				         PAR=PARmax-((PARmax-PARmin)*((double)iter2)/((double)Maxiter2)); 
						 if((iter2%2)==1) 
						 {
						 Xnew2[k]=MM2[L][k]+(((Math.random()-0.5)*2.0)*((Upper2[k]-Lower2[k])/200.0));	 
						 } 
						 else
						 {
						 Xnew2[k]=MM2[L][H]+(((Math.random()-0.5)*2.0)*((Upper2[k]-Lower2[k])/200.0)); 
						 }  
						   if(Math.random()<PAR)// pitch adjustment
						   {
							  for(int kk=0;kk<D;kk++)
						      {Xnew2[kk]=best3[kk];}	   
							   
						   } 
						   
				       }
					   else // Randomization
					   {
					      Xnew2[k]=Lower2[k]+(Math.random()*(Upper2[k]-Lower2[k]));	    
					   }
					   
					}	
			
			if(ff.func(Xnew2)<ff.func(worst3))		
			{
			    for(int i=0;i<D;i++)
			    {worst3[i]=Xnew2[i];}	
				
			}	
				
			for(int i=0;i<PMN;i++)
		    {fitness3[i]=ff.func(MM2[i]);}
		          
		    // get minimum index and best solution
		    double[] dd9=getminval_index(fitness3);
		    min2=dd9[0];
		    mindex2=(int)dd9[1];
		    for(int i=0;i<D;i++)
		    {best3[i]=MM2[mindex2][i];}
		          
		    // get maximum index and replace new worst value
		    double[] dd10=getmaxval_index(fitness3);
		    max2=dd10[0];
		    maxdex2=(int)dd10[1];
		    for(int i=0;i<D;i++)
		    {MM2[maxdex2][i]=worst3[i];}	
				
			 // Procedure for determining lower and upper bounds for every column
			  min=0.0;
			 for(int k=0;k<D;k++)
			 {
				 min=MM2[0][k];
			     for(int i=0;i<PMN;i++)
			     {
				  if(MM2[i][k]<min)
				  {min=MM2[i][k];}    
				 }
				 Lower2[k]=min;	 
			 }
			
			 
			 
			  max=0.0;
			 for(int k=0;k<D;k++)
			 {
				 max=MM2[0][k];
			     for(int i=0;i<PMN;i++)
			     {
				  if(MM2[i][k]>max)
				  {max=MM2[i][k];}    
				 }
				 Upper2[k]=max;	 
			 }	
				
			 iter2++;	
				
			}
			double[][] dp=new double[2][D];
			dp[0][0]=min2;
			for(int i=0;i<D;i++)
			{dp[1][i]=best3[i];}
			
			return dp;
			
		}
		
		double[] getmaxval_index(double a[])
	    {
	    double m=0.0;
	    double b[]=new double[a.length];
	    for(int i=0;i<a.length;i++)
	    {b[i]=a[i];}
	    double maxval=a[0];
	    for(int j=0;j<a.length;j++)
	    {if(a[j]>maxval){maxval=a[j];}}
	    for(int i=0;i<b.length;i++)
	    {if(b[i]==maxval){m=i;break;}}
	    double dep2[]=new double[2];
	    dep2[0]=maxval;
	    dep2[1]=m;
	    return dep2;	 
	    }
	    
	    
	    
	    double[] getminval_index(double[] a)
	    {
	    double m=0.0;	 
	    double b[]=new double[a.length];
	    for(int i=0;i<a.length;i++)
	    {b[i]=a[i];}  
	    double minval=a[0];
	    for(int i=0;i<a.length;i++)
	    {if(a[i]<minval){minval=a[i];}}
	    for(int i=0;i<a.length;i++)
	    {if(b[i]==minval){m=i;break;}};
	    double[] dep=new double[2];
	    dep[0]=minval;
	    dep[1]=m;
	    return dep;	 
	    }
	    
	    void toStringnew()
	    {
		    double out[][]=solution();
		    System.out.println("Optimized value = "+out[0][0]);
		    for(int i=0;i<D;i++)
		    {System.out.println("x["+i+"] = "+out[1][i]);}
		    
		}
	    
}