python 类鸟群Boids

import sys, argparse
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.spatial.distance import squareform, pdist, cdist  # 计算点之间的距离
from numpy.linalg import norm
from matplotlib.colors import ListedColormap
import matplotlib.patches as patches

width, height = 800, 800  # 设置屏幕上模拟窗口的宽度和高度


class Boids:
    def __init__(self, N):
        """ initialize the Boid simulation"""
        # 初始化位置和速度
        self.pos = [width / 2.0, height / 2.0] + 10 * np.random.rand(2 * N).reshape(N,
                                                                                    2)  # '''创建一个 numpy 数组 pos，对窗口中心加上 10 个单位以内的随机偏移。代码np.random.rand（2 * N）创建了一个一维数组，包含范围在[0，1]的 2N 个随机数。然后 reshape()调用将它转换成二维数组的形状（N，2），它将用于保存类鸟群个体的位置。'''
        # normalized random velocities
        angles = 2 * math.pi * np.random.rand(N)
        self.vel = np.array(list(zip(np.sin(angles), np.cos(angles))))
        self.N = N  # 生成一个数组，包含 N 个随机角度，范围在[0, 2pi]，
        # min dist of approach
        self.minDist = 25.0
        # max magnitude of velocities calculated by "rules"
        self.maxRuleVel = 0.03
        # max maginitude of final velocity
        self.maxVel = 2.0

    def tick(self, frameNum, pts, beak):
        """Update the simulation by one time step."""
        pdis_v=pdist(self.pos)
        self.distMatrix = squareform(pdist(self.pos))  # 用 squareform()和 pdist()方法来计算一组点之间两两的距离
        # apply rules:
        self.vel += self.applyRules()
        self.limit(self.vel, self.maxVel)
        self.pos += self.vel
        self.applyBC()
        # update data
        pts.set_data(self.pos.reshape(2 * self.N)[::2],
                     self.pos.reshape(2 * self.N)[1::2])
        vec = self.pos + 10 * self.vel / self.maxVel
        beak.set_data(vec.reshape(2 * self.N)[::2],
                      vec.reshape(2 * self.N)[1::2])

    def limitVec(self, vec, maxVal):
        """limit magnitide of 2D vector"""
        mag = norm(vec)
        if mag > maxVal:
            vec[0], vec[1] = vec[0] * maxVal / mag, vec[1] * maxVal / mag

    def limit(self, X, maxVal):
        """limit magnitide of 2D vectors in array X to maxValue"""
        for vec in X:
            self.limitVec(vec, maxVal)

    def applyBC(self):
        """apply boundary conditions"""
        deltaR = 2.0  # 该行中的deltaR提供了一个微小的缓冲区，它允许类鸟群个体开始从相反方向回来之前移出小块之外一点，从而产生更好的视觉效果
        for index, coord in enumerate(self.pos):
            if coord[0] > width + deltaR:
                coord[0] = - deltaR
            if coord[0] < - deltaR:
                coord[0] = width + deltaR
            if coord[1] > height + deltaR:
                coord[1] = - deltaR
            if coord[1] < - deltaR:
                coord[1] = height + deltaR
            """
            if ((coord[1] - 225) ** 2 + (coord[0] - 225) ** 2) <= 3600:
                self.vel[index] = -self.vel[index] * 2
            """

    def applyRules(self):
        # apply rule #1 - Separation
        D = self.distMatrix < 25.0
        vel=np.random.rand(self.N*2).reshape(self.N,2)
        vel = self.pos * D.sum(axis=1).reshape(self.N, 1) - D.dot(self.pos)
        self.limit(vel, self.maxRuleVel)

        # different distance threshold
        D = self.distMatrix < 50.0

        # apply rule #2 - 列队
        vel2 = D.dot(self.vel) # 50个单位内个体的速度向量和
        self.limit(vel2, self.maxRuleVel)
        vel += vel2;

        # apply rule #1 - 聚集
        vel3 = D.dot(self.pos) - self.pos
        self.limit(vel3, self.maxRuleVel)
        vel += vel3

        return vel

    def buttonPress(self, event):
        """event handler for matplotlib button presses"""
        # left click - add a boid
        if event.button == 1:
            self.pos = np.concatenate((self.pos,
                                       np.array([[event.xdata, event.ydata]])),
                                      axis=0)
            # random velocity
            angles = 2 * math.pi * np.random.rand(1)
            v = np.array(list(zip(np.sin(angles), np.cos(angles))))
            self.vel = np.concatenate((self.vel, v), axis=0)  # 拼接
            self.N += 1
            # right click - scatter
        elif event.button == 3:
            # add scattering velocity
            self.vel += 0.1 * (self.pos - np.array([[event.xdata, event.ydata]]))


def tick(frameNum, pts, beak, boids):
    # print frameNum
    """update function for animation"""
    boids.tick(frameNum, pts, beak)
    return pts, beak


# main() function
def main():
    # use sys.argv if needed
    print('starting boids...')

    parser = argparse.ArgumentParser(description="Implementing Craig Reynold's Boids...")
    # add arguments
    parser.add_argument('--num-boids', dest='N', required=False)
    args = parser.parse_args()

    # number of boids
    N = 100
    if args.N:
        N = int(args.N)

    # create boids
    boids = Boids(N)

    # setup plot
    fig = plt.figure(facecolor='pink')
    ax = plt.axes(xlim=(0, width), ylim=(0, height), facecolor='lightskyblue')
   # ax.add_patch(patches.Rectangle((200, 200), 50, 50, linewidth=1, edgecolor='b', facecolor='b'))
    pts, = ax.plot([], [], markersize=10,
                   c='k', marker='o', ls='None')
    beak, = ax.plot([], [], markersize=4,
                    c='r', marker='o', ls='None')

    anim = animation.FuncAnimation(fig, tick, fargs=(pts, beak, boids),
                                   interval=50)

    # add a "button press" event handler
    cid = fig.canvas.mpl_connect('button_press_event', boids.buttonPress)

    plt.show()

def test():
    pos=np.array([[1,2],[8,9],[2,3]])
    D=np.array([[1,0,1],[0,1,0],[1,0,1]])
    D_dot_pos= D.dot(pos)
    D1=np.array([[1,0,0],[1,0,0],[1,0,1]])
    D1_sum= D1.sum(axis=0);
    pass
# call main
if __name__ == '__main__':
    main()
    #test()

pos 是N*2 的位置向量

 self.distMatrix = squareform(pdist(self.pos)) 

pdist(self.pos) 计算向量两两之间的距离，返回数组，元素个数是 (N*N - N )/2  ， 距离矩阵是对称矩阵，对角线元素是0，pdis只保存上三角部分（不包括对角线）

squareform 将pdist返回的结果转换成N*N矩阵 D,里面的元素M(i,j) 表示i->j的距离， 由于对称 M(j,i)=M(i,j) 

--------------------------------------------

 D = self.distMatrix < 25.0 

返回距离值中小于25的1,0值矩阵，也是 N*N

D.dot(pos),  表示N*N 矩阵 乘以 N*2  结果也是N*2, 考虑列向量方式计算， D.dot(pos) 的意义是当前向量25范围内的向量之和。如当前范围内有 v1,v6,v8,那么这个结果是

v=[x1+x6+x8,y1+y6+y8]

 -------------------------------------

# apply rule #1 - Separation
D = self.distMatrix < 25.0
vel = self.pos * D.sum(axis=1).reshape(self.N, 1) - D.dot(self.pos)
self.limit(vel, self.maxRuleVel)

规则1：25个单位内 一小群个体间保持一定距离

D.sum(axis=1),是按行方向累加D中的行向量（由于对称，按列累加也一样），累加后转N*1矩阵，效果就是pos 对应向量  乘以累加的数， 比方一个小群内有6个鸟，那么累加完就是

pos*6,相当于6个这样的向量， 减去D.dot(pos)  （参考上面），相当于对应的位置向量分别减去小群内其他个体的位置向量的差之和，通过简单的画图可以发现这个向量迫使群内每个个体朝外（远离群）

----------------------------------

# apply rule #2 - 列队
vel2 = D.dot(self.vel)
self.limit(vel2, self.maxRuleVel)
vel += vel2;

这个规则以个体的速度为处理对象，即小群体内的速度向量之和， limit的作用是控制速度变化不能太过剧烈

---------------------------------

# apply rule #1 - 聚集
vel3 = D.dot(self.pos) - self.pos
self.limit(vel3, self.maxRuleVel)
vel += vel3

这个规则是当前小群体内所有位置向量之和 减掉[自身], 即指向合成向量的方向，效果是靠近小群体， 考虑pos向量时时更新，即合成向量时时更新，那么这个向量使个体靠近群体