//////////////////////////////////////////////////////////////////////
//                                                                  //
//   Cuckoo Search (CS) algorithm by Xin-She Yang and Suash Deb     // 
//    Programmer : Oğuz Emrah Turgut                                // 
//////////////////////////////////////////////////////////////////////
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.Random;
import java.util.ArrayList;

public class cuckoo_search
{	
	int nn=50;//number of nests
	double pa=0.35; // discovery_rate_of_alien_eggs
	double tol=1.0e-4;//tolerance
	int nd=2;
	
	double[] Lb;
	double[] Ub;
	double Lbbb=-5.0;
	double Ubbb=5.0;
		
	public double func(double x[]) 
    {
    //çözümü istenen fonksiyon	
     double ff=0;
    int n=x.length;
    for(int i=0;i<n;i++)
    {ff+=-Math.sin(x[i])*Math.pow(Math.sin((i+1)*x[i]*x[i]/Math.PI),20);}
     return ff; //maksimum testi
    }	
    public double getminval(double array[])
    {
	    double minval=array[0];    	    
	    for(int i=0;i<array.length;i++)
	    {
		    if(array[i]<minval)
		    {minval=array[i];}    		  
		}
	  return minval;	    
	} 
    
	
	public double gamma(double xx)
     {   //Gamma fonksiyonu
	double a[]=lngamma(xx);
 	return a[1]*Math.exp(a[0]);
     }   
	public  double[] lngamma(double x)
    {//logaritmik hata fonksiyonu
    //xx[0] hata fonksiyonunun logaritmasının değerini 
    //xx[1] hata fonksiyonunun kendisinin işaretini verir
	double xx[]=new double[2];    
    double result;
    double a;
    double b;
    double c;
    double p;
    double q;
    double u;
    double w;
    double z;
    int i;
    double logpi;
    double ls2pi;
    double tmp;
    double sgngam = 1;
    logpi = 1.14472988584940017414;
    ls2pi = 0.91893853320467274178;
    if( x<-34.0 )
    {
        q = -x;
        double yy[]=lngamma(q);
        w = yy[0];
        tmp=yy[1];
        p = Math.floor(q);
        i = (int)Math.round(p);
        if( i%2==0 )
        {
            sgngam = -1;
        }
        else
        {
            sgngam = 1;
        }
        z = q-p;
        if( z>0.5 )
        {
            p = p+1;
            z = p-q;
        }
        z = q*Math.sin(Math.PI*z);
        result = logpi-Math.log(z)-w;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    if( x<13 )
    {
        z = 1;
        p = 0;
        u = x;
        while(u>=3)
        {
            p = p-1;
            u = x+p;
            z = z*u;
        }
        while(u<2)
        {
            z = z/u;
            p = p+1;
            u = x+p;
        }
        if( z<0 )
        {
            sgngam = -1;
            z = -z;
        }
        else
        {
            sgngam = 1;
        }
        if( u==2 )
        {
            result = Math.log(z);
            //return result;
            xx[0]=result;
            xx[1]=sgngam;
            return xx;
        }
        p = p-2;
        x = x+p;
        b = -1378.25152569120859100;
        b = -38801.6315134637840924+x*b;
        b = -331612.992738871184744+x*b;
        b = -1162370.97492762307383+x*b;
        b = -1721737.00820839662146+x*b;
        b = -853555.664245765465627+x*b;
        c = 1;
        c = -351.815701436523470549+x*c;
        c = -17064.2106651881159223+x*c;
        c = -220528.590553854454839+x*c;
        c = -1139334.44367982507207+x*c;
        c = -2532523.07177582951285+x*c;
        c = -2018891.41433532773231+x*c;
        p = x*b/c;
        result = Math.log(z)+p;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    q = (x-0.5)*Math.log(x)-x+ls2pi;
    if( x>100000000 )
    {
        result = q;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    p = 1/(x*x);
    if( x>=1000.0 )
    {
        q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
    }
    else
    {
        a = 8.11614167470508450300*0.0001;
        a = -5.95061904284301438324*0.0001+p*a;
        a = 7.93650340457716943945*0.0001+p*a;
        a = -2.77777777730099687205*0.001+p*a;
        a = 8.33333333333331927722*0.01+p*a;
        q = q+a/x;
    }
    xx[0]= q;
    xx[1]=sgngam;
    return xx;
    }
	
    public double[][][] init() //initialize
	{
	 double[][]	nest=new double[nn][nd];
	 double[]	fitness=new double[nn];
		 Lb=new double[nd];//lower bounds
		 Ub=new double[nd];//upper bounds
		double[] randmat=new double[nd];
		double[] nestdepo=new double[nd];
		//applying upper and lower bounds
		for(int i=0;i<nd;i++)
		{
		 	Lb[i]=-5;
		 	Ub[i]=5;
		}
		// random initial solutions
	  for(int i=0;i<nn;i++)
	  { 
		for(int j=0;j<nd;j++)
		{randmat[j]=Math.random();
		nestdepo[j]=Lb[j]+((Ub[j]-Lb[j])*randmat[j]);
		nest[i][j]=nestdepo[j];
		}
      }	
	  
	  
	  for(int i=0;i<nn;i++)
	  {
	  fitness[i]=Math.pow(10,7); 	  
	  }
	 
	  double[][][] dep=new double[2][nn][nd];
	  
	  for(int i=0;i<nn;i++)
	  {
	      for(int j=0;j<nd;j++)
	      {
	           dep[0][i][j]=nest[i][j]; 	      
		  }	  
	  }
	  
	  for(int i=0;i<nn;i++)
	  {
	  dep[1][i][0]=fitness[i];	  
	  } 
	
	  return dep;
	}
	
	public double[][] get_cuckoos(double[][] nest,double[] best,double Lbb,double Ubb) // get cuckoos by random walk
	{   
		// Get cuckoos by random walk
		double[] s=new double[nd];
		double[] u=new double[nd];
		double[] v=new double[nd];
		double[] step=new double[nd];
		double[] stepsize=new double[nd];
		// introducing Levy flight
		int n1=nest.length;
		
		// Levy exponent and coefficient
		//For details, see equation (2.21), Page 16 (chapter 2) of the book
		//X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press,(2010)
		
		double beta=1.5;
		double sigma= Math.pow( (  gamma(1+beta)  *  Math.sin(3.1415*beta/2) / (gamma( (1+beta)/2 )*beta*Math.pow(2,((beta-1)/2))    )     ),(1/beta) );
		
		
		
		Random rndm=new Random();
		// to imply Uniformly distrubuted random numbers
		// Gaussian distribution
		
	  for(int j=0;j<nn;j++)	
	     {
	          s=nest[j];
		    for(int i=0;i<nd;i++)
		   {
			u[i]=(rndm.nextGaussian())*sigma;
			v[i]=(rndm.nextGaussian());
			step[i]=u[i]/Math.pow((Math.abs(v[i])),(1/beta) );
			stepsize[i]=0.01*step[i]*(s[i]-best[i]);
			s[i]=s[i]+stepsize[i]*(rndm.nextGaussian());
		   }
		      nest[j]=s;
         }   
		      double[][] exit=simplebounds(nest,Lbbb,Ubbb);  
          
         return exit;
		}
	
	public double[][] simplebounds(double[][] s,double Lb, double Ub)
	{
		// Application of simple constraints
		double[][] ns_temp=new double[s.length][s[0].length];
		
		for(int i=0;i<s.length;i++)
		{
		   for(int j=0;j<s[0].length;j++)
		   {
		      if(s[i][j]<Lb)
		      {s[i][j]=Lb;}
		       else if(s[i][j]>Ub)
		      {s[i][j]=Ub;}	   
		   }	
		}
		for(int i=0;i<s.length;i++)
		{
		   for(int j=0;j<s[0].length;j++)
		   {
		    ns_temp[i][j]=s[i][j];	   
		   }	
		}
	  
	    return ns_temp;
	  }
	
	// Find current best nest
	public double[][][] get_best_nest(double nest[][],double newnest[][],double[] fitness)
	{    
		double fnew;
		
		double []   best=new double[nd];
		for(int i=0;i<nn;i++)
		{
		   for(int j=0;j<nd;j++)
		   {
		      newnest[i][j]=nest[i][j];	   
		   }	
		}
		
		  // Evaluating all new solutions
		
	    for(int j=0;j<nn;j++)
	    {	  
		  for(int i=0;i<nd;i++)
		  {  
			 fnew=func(newnest[j]);
	       
		     if(fnew<fitness[j])
		     {
			  fitness[j]=fnew;
		      nest[j][i]=newnest[j][i];
		     }
		  }
	    }  
		  double min=getminval(fitness);
		  
		  int index=0;
		  for(int i=0;i<fitness.length;i++)
		   {
			 if(fitness[i]==min)
		 	      {break;} 
		     index++;
		   }
		  
		  
		  for(int i=0;i<nd;i++)
		   {
		   best[i]=nest[index][i];	  
		   }
		  
	
		  double[][][] depo =new double[4][nn][nd];
		  
		  for(int i=0;i<nn;i++){
		     for(int j=0;j<nd;j++){
		               depo[0][i][j]=nest[i][j];
	               }   
          }
	               
	      depo[1][0][0]=min;
	      
	      for(int i=0;i<nn;i++)
	      {depo[2][i][0]=fitness[i];}
          
	      for(int i=0;i<nd;i++)
	      {depo[3][1][i]=best[i];}
          return depo;
          }
	
	    
	
	    public double[] randperm(double[] firstarray)
	    { 
		    // random permutation of an array 
		    int n=firstarray.length;
		    Double[] db=new Double[n];
		    for(int i=0;i<n;i++)
		    {
			db[i]=firstarray[i];    
			}
		    List<Double> list = Arrays.asList(db);
		   
		    Collections.shuffle(list);
		    Object array[]=list.toArray();
		    double[] arrays=new double[array.length];
		    String str;
		    
		    for(int i=0;i<array.length;i++)
		    {
			  str=array[i].toString();  
			  arrays[i]=Double.valueOf(str).doubleValue();  
			} 
			
			return arrays;
			
		 } 
	     
	     
	     
	     // Replace some nest by constructing new solutions
	    public double[][] empty_nests(double[][] nest)
	    { 
		    // A fraction of worse nests are discovered with a probability pa
		   int n=nest.length; 
		   double[][] nest2=new double[nest.length][nest[0].length];
		   double[] temp11=new double[nest.length];
		   double[] temp22=new double[nest.length];
		   double[][] newnest=new double[nest.length][nest[0].length];
		   double[][] temp=new double[nest.length][nest[0].length];
		   double[][] temp3=new double[nest.length][nest[0].length];
		   double[][] temp4=new double[nest.length][nest[0].length];
		   double[] temp1=new double[nest.length];
		   double[] temp2=new double[nest.length];
		   double[][] K=new double[nest.length][nest[0].length];
		   
		   for(int i=0;i<nn;i++)
		   {
		      for(int j=0;j<nd;j++)
		      {
		         nest2[i][j]=nest[i][j];	      
			  }	   
		   }
		   
		    for(int i=0;i<nn;i++)
		   {
		      for(int j=0;j<nd;j++)
		      {
		         temp[i][j]=Math.random();	      
			        K[i][j]=0.0;
		      }
	       } 
	       
	       for(int i=0;i<nn;i++)
	       {
		      for(int j=0;j<nd;j++)
		      {
		        if( temp[i][j]<pa)
		        {K[i][j]=0.0;}
			    else if(temp[i][j]>pa)
			    {K[i][j]=1.0;}  
			  }    
		   }
	       
		for(int j=0;j<nd;j++)
		{   
	       for(int i=0;i<nn;i++)
	       {
		   temp1[i]=nest[i][j];    
		   }
		   temp2=randperm(temp1);
		   for(int i=0;i<nn;i++)
		   {
		    nest[i][j]=temp2[i]; 	   
		   }
	    }   
	    
	    for(int j=0;j<nd;j++)
		{   
	       for(int i=0;i<nn;i++)
	       {
		   temp11[i]=nest2[i][j];    
		   }
		   temp22=randperm(temp11);
		   for(int i=0;i<nn;i++)
		   {
		    nest2[i][j]=temp22[i]; 	   
		   }
	    }   
	    
	    double[][] subs=Matrix.substract(nest,nest2);
	    
		for(int i=0;i<nn;i++)
		{
		  for(int j=0;j<nd;j++)
		  {
			temp3[i][j]=Math.random()*subs[i][j];
			temp4[i][j]=temp3[i][j]*K[i][j];
		    newnest[i][j]=nest[i][j]+temp4[i][j];
			
		  }	
	    }   
	    return newnest;
	    
	    }
	     
	     
	public static void main(String args[])
	{
	/////////////////////////////////////////////////////////////	
	int nn=50;//number of nests
	double pa=0.35; // discovery_rate_of_alien_eggs
	double tol=1.0e-4;//tolerance
	int nd=2;//number of variables
	double Lbb=-5.0;
	double Ubb=5.0;
	double[][] nest=new double[nn][nd];
	double[] fitness=new double[nn];
	double[] best=new double[nd];
	double[] bestnest=new double[nd];
	double[][] newnest=new double[nn][nd]; 
	double[][][] d2;
	double fmin,fnew;
	int iter;
	cuckoo_search cs=new cuckoo_search();
	
	// initializing values 
	double[][][] d1=cs.init();
	for(int i=0;i<nn;i++)
	{      fitness[i]=d1[1][i][0];
	   for(int j=0;j<nd;j++)
	   {   
	       nest[i][j]=d1[0][i][j]; 	   
	   }	
    }
	
    
    
	
	d2=cs.get_best_nest(nest,nest,fitness);
	
	for(int i=0;i<nn;i++)
	{   
		fitness[i]=d2[2][i][0];
	    for(int j=0;j<nd;j++)
	    {     bestnest[j]=d2[3][1][j];
	          nest[i][j]=d2[0][i][j];	    
		}	
    }
	fmin=d2[1][0][0];
       
	iter=0;
	
	while(fmin>tol)
	{
		newnest=cs.get_cuckoos(nest,bestnest,Lbb,Ubb);
		
		d2=cs.get_best_nest(nest,newnest,fitness);
	
	for(int i=0;i<nn;i++)
	{   
		fitness[i]=d2[2][i][0];
	    for(int j=0;j<nd;j++)
	    {     best[j]=d2[3][1][j];
	          nest[i][j]=d2[0][i][j];	    
		}	
    }
	fnew=d2[1][0][0];
		
	iter+=nn;
	
	 newnest=cs.empty_nests(nest);
		
	d2=cs.get_best_nest(nest,newnest,fitness);
	
	for(int i=0;i<nn;i++)
	{   
		fitness[i]=d2[2][i][0];
	    for(int j=0;j<nd;j++)
	    {     best[j]=d2[3][1][j];
	          nest[i][j]=d2[0][i][j];	    
		}	
    }
	fnew=d2[1][0][0]; 
		
	iter+=nn; 
	
	if(fnew<fmin)
	{
	fmin=fnew;	
    bestnest=best;
	}	
	}
	System.out.println(Matrix.toString(bestnest));
	System.out.println(fmin);
	System.out.println(iter+" iterasyon");
    }
	
	
	
	
	
	
	
	
}