4.Julia数组、向量化和广播

数组的基本概念和用法不在这里多说了。对于科学计算，数组是不可或缺的。我们先写一个函数：

function f(x,y)return x*yend

再调用这个函数：

print(f(3,5))

运行这个函数得：

现在我们有A,B两数组，元素一样多。我们想把这两个数组逐个元素相乘，可以写一个for循环来完成，也可以利用向量化的方法来完成。向量化就是函数作用在数组的每个元素上。代码如下：

A=range(1,9,9) #[1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0]B=range(1,9,9) #[1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0]C=f.(A,B)print("$C\n")#运行结果：[1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0]

用向量化的方式替代for循环的好处是代码简洁，含义明确。即函数逐个作用在元素上，且要求两个数组大小一样。

如果想A、B中每个元素滚动相乘。像乘法口诀表一样，这需要两个嵌套的for循环：

print( "$A\n$B\n")for i=1 : length(A)for j =1 : length(B)C[i,j]=f(A[i],B[j])endendprint(C)println()#结果如下
[1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0;2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0;3.0 6.0 9.0 12.0 15.0 18.0 21.0 24.0 27.0;4.0 8.0 12.0 16.0 20.0 24.0 28.0 32.0 36.0;5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0;6.0 12.0 18.0 24.0 30.0 36.0 42.0 48.0 54.0;7.0 14.0 21.0 28.0 35.0 42.0 49.0 56.0 63.0;8.0 16.0 24.0 32.0 40.0 48.0 56.0 64.0 72.0;9.0 18.0 27.0 36.0 45.0 54.0 63.0 72.0 81.0]

那可不可以像向量化一样，替代嵌套的for循环呢？答案是利用广播。广播的作用可以理解为不同维度的数组之间的运算。

[1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0]

先把数组A进行转置：

D=reshape(A,(1,:))

print(f.(D,B))

结果正是我们想要的。

好了，回到我的问题：一个飞机场每天可能会有上万条飞机飞行轨迹，每条飞行轨迹有约30个点组成，约十万个计算点。

这正是三个for循环嵌套。可以利用三个不同维度的数组广播来实现。

# 空间坐标
struct Point3Dx::Float64y::Float64z::Float64 # aeroplane height
end# 飞行剖面
struct ProfilePoint::Point3D  # 空间位置V::Float64  # groundspeed。速度s::Float64  # distance along the ground track。点到起点的距离P::Float64  # noise-related power parameter。功率ϵ::Float64 # bank angle。倾斜角end
# 飞行轨迹段
struct SegmentProfilebegainPf  #段的起点参数endPf  #段的终点参数
end function getsegment(pf::AbstractArray) SP=[]for i =1 : length(pf)-1insert!(SP,i,SegmentProfile(pf[i],pf[i+1]))endreturn SP
end
# 计算空间两点的距离
function disdance(point1::Point3D,point2::Point3D)::Float64return sqrt((point2.x-point1.x)^2+(point2.y-point1.y)^2+(point2.z-point1.z)^2)
end
# 点到线的垂点
function GetFootOfPerpendicular(recepter::Point3D,Point1::Point3D,point2::Point3D)dx = Point1.x - point2.xdy = Point1.y - point2.ydz = Point1.z - point2.zif abs(dx) < 0.00000001 && abs(dy) < 0.00000001 && abs(dz) < 0.00000001 return Point1end u = (recepter.x - Point1.x)*(Point1.x - point2.x) +(recepter.y - Point1.y)*(Point1.y - point2.y) + (recepter.z - Point1.z)*(Point1.z - point2.z)u = u/((dx*dx)+(dy*dy)+(dz*dz))return Point3D(Point1.x + u*dx, Point1.y + u*dy, Point1.z + u*dz)end# 计算飞行剖面段和观察者之间的参数Figure 4-2
function getgeometic(recepter::Point3D,SP::SegmentProfile)#Recept_Sp=GetFootOfPerpendicular(recepter,SP.begainPf.Point,SP.endPf.Point) # 是线段上垂直于最接近观察者的点或的延伸点#= d1=disdance(SP.begainPf.Point,recepter)# 段开始和观察者之间的距离d2=disdance(SP.endPf.Point,recepter) # 段结束和观察者之间的距离q=disdance(SP.begainPf.Point,Recept_Sp)  # 是S1到Sp的距离（如果观察者位置在段后面，则为负）λ=disdance(SP.begainPf.Point,SP.endPf.Point)  # 是飞行路径段的长度if disdance(SP.endPf.Point,Recept_Sp)>disdance(SP.begainPf.Point,Recept_Sp) &&  disdance(SP.endPf.Point,Recept_Sp)>λq=-abs(q)endds=0 # 是观察者和线段之间的最短距离dp=disdance(recepter,Recept_Sp) # 是观察者与扩展段之间的垂直距离（最小倾斜范围）P=0P1=SP.begainPf.PP2=SP.endPf.Pif q<0ds=d1P=P1elseif 0<=q && q<=λds=dpP=sqrt(P1^2+(q/λ)*(P2^2-P1^2))elseds=d2P=P2endreturn (dp,ds,P) =#return 0 end#测试代码如下：
# 构造飞行剖面数据
listpf=[Profile(Point3D(0,0,0),0,0,0,0),Profile(Point3D(300,0,0),100,300,50,0),Profile(Point3D(600,0,0),130,600,100,0),Profile(Point3D(1000,0,0),150,1000,150,0),Profile(Point3D(1500,0,0),200,1500,200,0)]#print(listpf)
#print("\n")
listSg=getsegment(listpf) # 根据飞行剖面数据获得飞行段数据
print(listSg)
print("\n")
Recept=[Point3D(400,-200,0),Point3D(700,-300,0),Point3D(1200,-400,0)]
for Sg in listSg #循环嵌套for r in ReceptRes=getgeometic(r,Sg)print(Res)print("\n")end
end
println()print(getsegment.(reshape(Recept,(1,:)),listSg)) #广播代替循环

运行结果错误！没能正确的实现广播。调试跟踪一下，没找到。

于是，我又得到一个原则，对于复杂的自定义结构数据，广播不一定好用，而且调试十分困难。这时候for循环不香吗？按我的需求10000*30*100000=300亿个数据，如果用广播会不会崩盘（这还不算每年有365天）？矢量化和广播还是用在它该用的地方吧。

我突然想起一本书：《流畅的Python》。Python那优美的代码让人心动。可惜我的脑子里有十几种语言搅在一起，有些地方思维已成定式。for循环在脑子里已经扎根，要我放弃他吗？先等等。如果Python的for循环非常快，它将是我的第一选择！