迭代函数系统 - iterated function system是一种创建分形图案的简单算法。下面我们用迭代函数系统来凭"空"生成一片树叶。下面的4组线性函数均可以根据一个二维平面点的点坐标(xi,yi)计算得到一个新的点坐标(xi+1,yi+1)。

知识产权协议

允许以教育/培训为目的向学生或受众进行免费引用，展示或者讲述，无须取得作者同意。

不允许以电子/纸质出版为目的进行摘抄或改编。

函数组0，选择概率0.01：

函数组1，选择概率0.07:

函数组2，选择概率0.07:

函数组3，选择概率0.85:

整个计算过程也是迭代的，我们首先选择坐标原点(0,0)作为(xi,yi)，选择上述函数组中的一个生成新坐标(xi+1,yi+1)，然后再把新坐标作为(xi,yi)，选择上述函数组中的一个计算得到下一个平面坐标点。在重复10万次后，我们就得到了平面上10万个点，将这10万个点在平面上画出来，即得该迭代函数系统的分形图案。

那么，在迭代时，选择哪个函数组进行迭代计算呢？答案是在符合概率要求的条件下随机选择。我们为每个函数组指定了选择概率，分别是1%，7%，7%和85%。可以看到，函数组3被选择的概率最高。

这次我们先看看绘图结果，如下。Amazing! 在只有数个参数的情况下，IFS成功地构造了树叶。读者可以点击左下方的放大镜，然后按住鼠标左键框选放大，可以看到，树叶中的某片子树叶，跟整片树叶一模一样！图像表现出高度的自相似性。

接下来分析代码。下面使用到的numpy是用C语言书写的多维数组模块，它是scipy等科学计算模块的基础。读者可以把numpy简单看作是多维的速度更快的列表(list)。numpy包可以在Windows命令行或者Linux终端中使用pip install numpy命令进行安装。

#IfsLeaf.pyimport numpy as npfrom matplotlib import pyplot as pltdef ifs(p,eq,n):    p = np.cumsum(p)              #对概率列表进行累加运算    rands = np.random.rand(n)          #生成n个值域为[0,1)的随机数    select = np.searchsorted(p,rands)      #生成方程组选择数组    result = np.zeros((n,2),dtype=np.float)    #n行2列的结果数组，元素值初始化为0    c = np.zeros(n,dtype=np.float)        #n个0值构成的c数组，用于表示点的颜色    #值为[0,0,1]的当前点坐标数据，可以视为3行1列的矩阵    pos = np.array([0,0,1],dtype=np.float)        for i in range(n):        eqidx = select[i]        tmp = np.dot(eq[eqidx],pos)        pos[:2] = tmp        result[i] = tmp        c[i] = eqidx    return result[:,0],result[:,1],c#方程组参数矩阵equations = [np.array([[0,0,0],[0,0.15,0]]),             np.array([[0.21, -0.19, 0], [0.24, 0.27, 1.5]]),             np.array([[-0.14, 0.26, 0], [0.26, 0.25, 0.47]]),             np.array([[0.87, 0, 0], [-0.05, 0.84, 1.54]])]x,y,c = ifs([0.01,0.07,0.07,0.85],equations,100000)fig,axes = plt.subplots(1,2,figsize=(6,5))     #生成一行2列两个子图，返回图及子图列表axes[0].scatter(x,y,s=1,c='g',marker='s',linewidths=0)  #以绿色画散点图axes[1].scatter(x,y,s=1,c=c,marker='s',linewidths=0)  #以不同颜色画散点图for ax in axes:    ax.axis("off")                    #不显示子图的坐标系plt.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)  #调整子图间距plt.show()                        #显示图表

我们使用np.array()函数转换生成了4个2行3列的矩阵，或者2维数组。其中，第0行用于计算xi+1，第1行用于计算yi+1。这4个矩阵依序存储于equations列表中。

函数ifs用于迭代计算，参数p为包含各方程组被选择概率的列表，参数eq对应代码中的equations列表，n则为指定的迭代次数 - 即需要迭代生成的平面坐标点的数量。

接下来我们得研究一下"奇怪"的select数组，它帮助我们按照约定概率随机选择迭代方程组。

p = np.cumsum(p)              #对概率列表进行累加运算 rands = np.random.rand(n)          #生成n个值域为[0,1)的随机数 select = np.searchsorted(p,rands)      #生成方程组选择数组

np.cumsum()对概率列表p进行累加运算。该函数执行前，p = [0.01, 0.07, 0.07, 0.85]，以p为参数执行完p = np.cumsum(p)以后，p值为array([0.01, 0.08, 0.15, 1. ])。显然，cumsum()函数就是逐一计算序列的前n项和。这个工作，我们在Python内使用一个for循环也可以完成，但是，Python是解释执行的“慢”语言，为了提高效率，我们使用了np.cumsum()函数。可以想象，这个由C语言写就的函数内部也会有一个循环，但由于它是C语言写的，所以执行效率要高得多!

np.random.rand(n)生成n个值域为[0,1)的随机数。为了便于说明问题，作者取n=8,生成了下述rands数组：

array([0.00165618, 0.32583924, 0.17917682, 0.94870972, 0.63444525,       0.14592743, 0.18000183, 0.16312187])

然后执行select = np.searchsorted(p,rands)，select的值如下：

array([0, 3, 3, 3, 3, 2, 3, 3], dtype=int64)

注意，此时的p值为array([0.01, 0.08, 0.15, 1. ])，它是一个递增的有序数组。searchsorted()函数逐一取出rands中呈均匀分布的随机数，比如0.00165618，然后在p数组中进行搜索，发现它小于p[0] = 0.01，导出结果下标0；下一个随机数为0.32583924，searchsorted()搜索发现其值大于p[2]，小于p[3]，导出结果下标3；依此类推，所有rands中的随机数均会产生一个搜索结果下标。显然，这些结果下标对应rands中的元素插入有序p数组中的“插入”位置。

上述select数组中的下标个数等于rands中值的个数，也就是n。select[i] = j意味着第i轮迭代将选择第j个方程组进行计算。简单分析可以看到，由于第3个方程组对应的概率累积值区间为[0.15,1)，因此，均匀分布的随机数落在这个区间的概率就是 1 - 0.15 = 0.85!

result = np.zeros((n,2),dtype=np.float)生成了一个n行2列的矩阵，其元素值的类型为np.float。同时，np.zeros()函数还把全部元素初始化为0。result二维数组存储了n个平面点坐标，第i个点的坐标在第i行。其中，result[i][0]是该点的x坐标，result[i][1]是该点的y坐标。

c = np.zeros(n,dtype=np.float)生成了一个由n个0值构成的c数组，c[i]表示第i个平面点的描绘颜色。

pos = np.array([0,0,1],dtype=np.float) 表示当前平面点的坐标，pos[0]为x坐标, pos[1]为y坐标, pos[2] = 1是根据方程计算的需要而添加的。pos数组有3个元素，可以视为3行1列的矩阵。

    for i in range(n):        eqidx = select[i]          #该次迭代用方程组编号        tmp = np.dot(eq[eqidx],pos)      #方程组系数矩形 x pos矩阵        pos[:2] = tmp            #更新当前平面点坐标，以备下轮迭代用        result[i] = tmp            #保存第i个平面点的坐标至结果数组        c[i] = eqidx            #方程组编号作为平面点颜色

接下来就是循环迭代计算。每次循环，先从select[i]取得当次迭代用的方程组编号eqidx，根据这一编号，eq[eqidx]返回对应编号的方程组的系数矩阵。如前所述，这个矩阵应为2行3列。同时，pos数组为3行1列的矩阵。根据线性代数常识，一个2行3列的矩阵与一个3行1列的矩阵可乘，np.dot()函数即执行该矩阵乘法。结果对象tmp应为一个2行1列的矩阵。读者可以自行验算一下，该矩阵乘法的执行过程就是把pos[0], pos[1]当作xi，yi并根据方程组计算xi+1，yi+1的过程。矩阵乘法涉及三重循环，由于np.dot()函数是用C语言书写的，其执行速度远远超过我们在Python中直接进行循环运算。

上述矩阵乘法结束后，tmp数据存储了新的平面点的坐标，tmp[0]为x坐标，tmp[1]为y坐标。pos[:2] = tmp的等号前边相当于给数组作了一个切片，这行代码等价于pos[0] = tmp[0], pos[1] = tmp[1]。这事实上更新了pos数据中的当前平面点坐标，同时保持pos[2] = 1不变，为下一轮迭代做好准备。result[i] = tmp将计算得到的第i个平面点坐标保存至result数组。c[i] = eqidx简单地将第i个平面点的对应方程组编号设置为该平面点描绘颜色。

ifs()函数返回了result[:,0],result[:,1],c三个1维数组，其长度均为n。从左到右，依次为n个平面点的x坐标, n个平面点的y坐标, n个平面点的颜色。result[:,0]可以视为对这个二维数组的切片，该切片参数表明选择全部0轴取值(即所有行），而1轴只取0值（即只取第0列）。

x,y,c = ifs([0.01,0.07,0.07,0.85],equations,100000)通过调用ifs函数迭代计算了10万个平面点。结果数据中，x[i]表示第i个平面点的x坐标，y[i]表示第i个平面点的y坐标，c[i]表示第i个平面点的描绘颜色。

最后就是数据展现部分。fig,axes = plt.subplots(1,2,figsize=(6,5)) 创建一个6英寸长，5英寸高的图(figure)，并在图中创建1行2列两个子图(axes)。该函数返回包含figure及子图列表的元组。axes[0].scatter(x,y,s=1,c='g',marker='s',linewidths=0)在0号子图上画散点图，颜色统一为绿色 - c='g'；点的形状为方形 - marker='s'；linewidths=0表明点的边线宽度为0，即不画边线；s = 1表明点的描绘尺寸。1号子图上的散点图与0号子图只有一个区别，就是每个平面点的颜色指定为该平面点生成时对应的方程组编号 - c = c。通过颜色，我们可以推断树叶中的哪些部分是由哪一个方程组生成出来的。

plt.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)调整了子图之间的间距等信息。

读者可以观察执行结果中的右图的颜色，可以看出：概率为1%的函数组0生成了叶杆部分(深红)、概率为7%的函数组1，2生成了左右两片子叶（蓝、绿）、概率为85%的函数组3则生成了上面所有的黄色部分。

需要说明的是，如果读者的matplotlib版本与作者测试时的版本不同，上述颜色可能会有差异。本小节的内容用代码实现很容易，但要真正理解为什么IFS的迭代方程组可以生成上述图像则相对困难，这需要更多的线性代数以及2维仿射变换 - affine transformation的知识，超出了本书的描述范围。

后记

本案例安排在教科书《Python编程基础及应用》的靠后部分，此时，读者已经具备了numpy的基础知识，故上述案例以numpy为基础构建。

事实上，上述案例不使用numpy， 直接用Python里的列表也是可以的。各位同行可以自己改改，如果不想改，可以发邮件问我要，微信也可以。

本案例节选自作者编写的教材及配套实验指导书。

《C++编程基础及应用》（高等教育出版社，出版过程中）

《Python编程基础及应用》，高等教育出版社

《Python编程基础及应用实验教程》，高等教育出版社

高校教师同行如果期望索取样书，教学支持资料，加群，请私信作者，联系时请提供学校及个人姓名为盼，各高校在读学生勿扰为谢。

青少年读者们如果期望系统性地学习Python及C/C++程序设计语言，欢迎尝试下述今日头条（西瓜）免费视频课程。

C/C++从入门到放弃（重庆大学现场版）

Python编程基础及应用（重庆大学现场版）