布谷鸟算法

import numpy as np
import scipy.special as sc_special
import matplotlib.pyplot as plt             #matplotlib的pyplot模块一般是最常用的，可以方便用户快速绘制二位图表
from matplotlib import cm 					 #matplotlib是python最著名的绘图库
from mpl_toolkits.mplot3d import Axes3D    #3D绘图

#绘制图像
def plot_3d(ax):
    x = np.arange(-3, 3, 0.1)       #在指定的间隔内返回均匀间隔的数字
    y = np.arange(-3, 3, 0.1)
    x,y = np.meshgrid(x, y)
    z = 3*(1-x)**2*np.e**(-x**2-(y+1)**2) - 10*(x/5-x**3-y**5)*np.e**(-x**2-y**2) - (np.e**(-(x+1)**2-y**2))/3
    ax.plot_surface(x,y,z,rstride=1,cstride=1,cmap=cm.coolwarm)# rstride:行之间的跨度  cstride:列之间的跨度,cmap是颜色映射表
    ax.set_zlim(-10,10)               ##坐标系的下边界和上边界
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    plt.pause(3)                      #延时
    plt.show()                        #显示

def fit_func(nest):
        x, y = nest
        return 3*(1-x)**2*np.e**(-x**2-(y+1)**2) - 10*(x/5-x**3-y**5)*np.e**(-x**2-y**2) - (np.e**(-(x+1)**2-y**2))/3

#计算适应度
def calc_fitness(fit_func, nests):
    n, m = nests.shape
    fitness = np.empty(n)

    for each_nest in range(n):
        fitness[each_nest] = fit_func(nests[each_nest])

    return fitness


"""
    Cuckoo search function
    ---------------------------------------------------
    Input parameters:
        n: 巢的数量
        m: 维数
        fit_func: 适用度函数
        lower_boundary: 搜索下边界
        upper_boundary: 搜索上边界
        iter_num: 迭代次数 (默认: 100)
        pa: 被宿主发现蛋的概率 (default: 0.25)
        beta:与问题规模相关的步长比例因子 (note: 1 < beta < 2) (default: 1.5)
        step_size:  与问题规模相关的步长比例因子 (default: 0.1)
    Output:
        最佳解决方案及值
"""

# 算法实现 
# cuckoo_search的目的是生成'最好巢穴位置'和'最优适应度'
def cuckoo_search(n, m, fit_func, lower_boundary, upper_boundary, iter_num = 100,pa = 0.25, beta = 1.5, step_size = 0.1):
    # 得到最初巢穴的位置
    nests = generate_nests(n, m, lower_boundary, upper_boundary) #生成巢穴位置
    fitness = calc_fitness(fit_func, nests)                      #计算适应度
    # 初始化巢的位置
    best_nest_index = np.argmax(fitness)
    best_fitness = fitness[best_nest_index]
    best_nest = nests[best_nest_index].copy()

    for _ in range(iter_num):
        #获取新的巢穴位置并用好的适应值替换掉不好的适应值
        nests = update_nests(fit_func, lower_boundary, upper_boundary, nests, best_nest, fitness, step_size)
        #在一定限制条件下蛋被丢弃寻找新巢
        nests = abandon_nests(nests, lower_boundary, upper_boundary, pa)    
        #计算适应度
        fitness = calc_fitness(fit_func, nests)

        x=np.empty((n, 1))
        y=np.empty((n, 1))
        z=np.empty((n, 1))
        for i in range(20):
           x[i],y[i]=nests[i]
        z=3*(1-x)**2*np.e**(-x**2-(y+1)**2) - 10*(x/5-x**3-y**5)*np.e**(-x**2-y**2) - (np.e**(-(x+1)**2-y**2))/3
        if 'sca' in locals():
            sca.remove()
        sca = ax.scatter(x, y, z, c='black', marker='o');
        plt.show();
        plt.pause(0.1)        #ax.scatter特征值散点图

        max_nest_index = np.argmax(fitness)
        max_fitness = fitness[max_nest_index]
        max_nest = nests[max_nest_index]

        if (max_fitness > best_fitness):
            best_nest = max_nest.copy()
            best_fitness = max_fitness

    return (best_nest, best_fitness) #返回最好的巢穴位置和最优适应度

'''
        n: 巢的数量
        m: 维数
        fit_func: 适用度函数
        lower_boundary: 下边界
'''
#生成巢穴位置
def generate_nests(n, m, lower_boundary, upper_boundary):
    lower_boundary = np.array(lower_boundary)  #转成矩阵形式
    upper_boundary = np.array(upper_boundary)
    nests = np.empty((n, m))                   #生成n行m列的数组

    for each_nest in range(n):
        nests[each_nest] = lower_boundary + np.array([np.random.rand() for _ in range(m)]) * (upper_boundary - lower_boundary) #生成在[-3,3]范围内的下x,y的n个样本，即20*2的矩阵

    return nests


#获取新的巢穴位置并用好的替换掉旧的不好的
def update_nests(fit_func, lower_boundary, upper_boundary, nests, best_nest, fitness, step_coefficient):
    lower_boundary = np.array(lower_boundary)
    upper_boundary = np.array(upper_boundary)
    n, m = nests.shape
    # 使用莱维飞行产生步长
    steps = levy_flight(n, m, 1.5)  #下面有定义该函数
    new_nests = nests.copy()

    for each_nest in range(n):
        # coefficient 0.01 is to avoid levy flights becoming too aggresive
        # and (nest[each_nest] - best_nest) could let the best nest be remained
        step_size = step_coefficient * steps[each_nest] * (nests[each_nest] - best_nest)
        step_direction = np.random.rand(m)
        new_nests[each_nest] += step_size * step_direction
        # apply boundary condtions
        new_nests[each_nest][new_nests[each_nest] < lower_boundary] = lower_boundary[new_nests[each_nest] < lower_boundary]
        new_nests[each_nest][new_nests[each_nest] > upper_boundary] = upper_boundary[new_nests[each_nest] > upper_boundary]

    new_fitness = calc_fitness(fit_func, new_nests)
    nests[new_fitness > fitness] = new_nests[new_fitness > fitness]

    return nests


#卵被丢弃寻找新巢
def abandon_nests(nests, lower_boundary, upper_boundary, pa):
    lower_boundary = np.array(lower_boundary)
    upper_boundary = np.array(upper_boundary)
    n, m = nests.shape
    for each_nest in range(n):
        if (np.random.rand() < pa):
            step_size = np.random.rand() * (nests[np.random.randint(0, n)] - nests[np.random.randint(0, n)])
            nests[each_nest] += step_size
            # apply boundary condtions
            nests[each_nest][nests[each_nest] < lower_boundary] = lower_boundary[nests[each_nest] < lower_boundary]
            nests[each_nest][nests[each_nest] > upper_boundary] = upper_boundary[nests[each_nest] > upper_boundary]

    return nests


#计算莱维飞行
def levy_flight(n, m, beta):
    sigma_u = (sc_special.gamma(1+beta)*np.sin(np.pi*beta/2)/(sc_special.gamma((1+beta)/2)*beta*(2**((beta-1)/2))))**(1/beta)
    sigma_v = 1

    u =  np.random.normal(0, sigma_u, (n, m))
    v = np.random.normal(0, sigma_v, (n, m))

    steps = u/((np.abs(v))**(1/beta))  #steps为20*2的矩阵

    return steps

#主函数
if __name__=='__main__':
    fig = plt.figure()    #绘制背景图
    ax = Axes3D(fig,auto_add_to_figure=False)
    fig.add_axes(ax)
    plt.ion()#将画图模式改为交互模式，程序遇到plt.show不会暂停，而是继续执行
    plot_3d(ax)

    best_nest, best_fitness = cuckoo_search(25, 2, fit_func, [-3, -3], [3, 3], step_size = 0.4)
    # 利用cuckoo_search生成'最好巢穴位置'和'最优适应度'
    # 展示出搜索到的最优结果坐标
    print('最大值为:%.5f, 在(%.5f, %.5f)处取到!'%(best_fitness, best_nest[0], best_nest[1]))
    plt.ioff()    #将画图交互模式关闭
    plot_3d(ax)