import numpy as np
import random
import math
import matplotlib.pyplot as plt


para=[20,1000,0.5,0.5]
fun_number=1   # 1 Griewank, 2 Rosenblock, 3 Easom 
n=para[0]      # Population size, typically 10 to 40
N_gen=para[1]  # Number of generations
A=para[2]      # Loudness  (constant or decreasing)
r=para[3]      # Pulse rate (constant or decreasing)
Qmin=0         # Frequency minimum
Qmax=2         # Frequency maximum
N_iter=0       # Total number of function evaluations
d=10           # Number of dimensions 

if(fun_number==1):
    optimum = 0
    xmin=-600
    xmax= 600
    def function(x):
        fr = 4000
        n=len(x)	
        s = 0
        p = 1
        for j in range(n): 
            s = s+x[j]*x[j]
        for j in range(n): 
            p = p*np.cos(x[j]/np.sqrt(j+1))
        return s/fr-p+1

elif(fun_number==2):
     optimum = 0
     xmin=-5
     xmax=10
     def function(x):
         n=len(x)
         s=0
         for j in range(n-1):
             s=s+100*(x[j]*x[j]-x[j+1])*(x[j]*x[j]-x[j+1])+(x[j]-1)*(x[j]-1)
         return s

elif(fun_number==3):
    optimum = -1
    xmin=-100
    xmax=100
    n=2
    def function(x):
        return -np.cos(x[0])*np.cos(x[1])*np.exp(-(x[0]-np.pi)*(x[0]-np.pi)-(x[1]-np.pi)*(x[1]-np.pi))


def LowerLimitBounds(d):
    Lb=-2*np.ones((1,d))
    print(Lb)
    return Lb


def UpperLimitBounds(d):
    Ub=-2*np.ones((1,d))
    print(Ub)
    return Ub


def Frequency(n):
    Q=np.zeros((n,1))
    print(Q)
    return Q


def Velocities(n,d):
    v=np.zeros((n,d))
    print(v)
    return v


def Solution(n):
    Sol=np.zeros((n,1))
    print(Sol)
    return Sol


def InitializeThePopulation(Lb,Ub):
    for i in range(n):
        Sol[i,:]=Lb+(Ub-Lb)*np.random.rand(1,d)
        Fitness[i]=Sol[i,:]
        return Fitness

def FindTheInitialBestPopulation(Fitness):
    [fmin,I]=np.min(Fitness)
    best=Sol[I,:]
    return best


def simplebounds(s,Lb,Ub):
  #Apply the lower bound vector
    ns_tmp=s
    I=ns_tmp<Lb
    ns_tmp[I]=Lb[I]
  #Apply the upper bound vector 
    J=ns_tmp>Ub;
    ns_tmp[J]=Ub[J]
  #Update this new move 
    s=ns_tmp
    return s

#Sphere function
def fun(u):
    z=np.sum(u**2)
    return z

Lb=LowerLimitBounds(d)
Ub=UpperLimitBounds(d)
Q=Frequency(n)
v=Velocities(n,d)
Sol=Solution(n)
Fitness=InitializeThePopulation(Lb,Ub)

"""
Bat Algorithm (essential part)
"""
# Start the iterations
for t in range(N_gen):
    #Loop over all bats/solutions
    for i in range(n):
        Q[i]=Qmin+(Qmin-Qmax)*np.random.rand()
        v[i,:]=v[i,:]+[Sol[i,:]-best]*Q[i]
        S[i,:]=Sol[i,:]+v[i,:]
        #Apply simple bounds/limits
        Sol[i,:]=simplebounds(Sol[i,:],Lb,Ub)
        if np.random.rand()>r:
            #The factor 0.001 limits the step sizes of random walks 
            S[i,:]=best+0.001*np.random.randint(1,d) 
        # Evaluate new solutions
        Fnew=Fun(S[i,:])
        # Update if the solution improves, or not too loud
        if (Fnew<=Fitness(i)) & (np.random.rand()<A):
            Sol[i,:]=S[i,:]
            Fitness[i]=Fnew
        # Update the current best solution
        if Fnew<=fmin:
            best=S[i,:]
            fmin=Fnew
    N_iter=N_iter+n
                
         
       
print("Number of evaluations: ",N_iter)
print("Best = ",best, "fmin = ",fmin)