4 复合数据类型进阶

这一部分介绍Julia的数组和数组运算, 以及收集容器(collections)类型的一般用法。

4.1 一维数组

Julia支持一维和多维的数组, 当一维数组的元素是数值时, 也可以理解成数学中的向量。

在程序中直接定义一个向量, 只要用方括号内写多个逗号分隔的数值,如

v1 = [2, 3, 5, 7, 11, 13, 17]
7-element Vector{Int64}:
  2
  3
  5
  7
 11
 13
 17
v2 = [1.5, 3, 4, 9.12];
@show v2;
## v2 = [1.5, 3.0, 4.0, 9.12]

其中v1是整数型的向量, v2是浮点型Float64的向量。 也可以定义元素为字符串的数组, 元素为不同类型的数组, 等等:

v3 = ["苹果", "桔子", "香蕉"];
@show v3;
## v3 = ["苹果", "桔子", "香蕉"]
v4 = [123, 3.14, "数学", [1, 2, 3]];
@show v4;
## v4 = Any[123, 3.14, "数学", [1, 2, 3]]

length(x)求向量x的元素个数,如

length(v1)
## 7

4.1.1 向量下标

4.1.1.1 访问单个元素

x是向量,i是正整数, x[i]表示向量的第i个元素。 第一个元素的下标为1,这种规定与R、FORTRAN语言相同, 但不同于Python、C、C++、JAVA语言。 如

v1 = [2, 3, 5, 7, 11, 13, 17]; v1[2]
## 3

end表示最后一个元素位置,如:

v1[end]
## 17
v1[end-1]
## 13

用begin表示第一个元素位置,如:

v1[begin]
## 2

为什么不直接用1表示第一个元素位置而是提供一个begin关键字? 这是因为Julia支持一种OffsetArray, 可以使得下标不必要从1开始计数。 在可能使用了OffsetArray的情况下, 可以在下标中用begin表示第一个元素的下标。

4.1.1.2 修改元素值

对元素赋值将在原地修改元素的值,如

v1[2] = 0
@show v1;
## v1 = [2, 0, 5, 7, 11, 13, 17]

这说明数组是“可变类型”(mutable), 即其中的元素可以原地修改。 字符串和元组则属于不可变类型(immutable)。

4.1.1.3 用范围作为下标

下标可以是一个范围,如

v1 = [2, 3, 5, 7, 11, 13, 17]
v1[2:4]
3-element Vector{Int64}:
 3
 5
 7

从python语言转移过来的读者要注意, Julia用1作为第一个下标, 而且下标2:4就表示第2,3,4号元素, 而在Python中2:4意思是\(2 \leq i < 4\), 代表第2, 3号元素(从0开始编号)。

在这种范围中,仍可用end表示最后一个下标, 用begin表示第一个下标,如:

v1[begin+2:end-2] |> show
## [5, 7, 11]

begin+2:end-2的范围是第3个到倒数第3个的范围。 从上面的例子也可以看出, :运算符的优先级要低于四则运算的优先级, 在其它编程语言中可能要写成(begin+2):(end-2), Julia可以写成begin+2:end-2, 但为可读性也可以写成(begin+2):(end-2)

范围可以是start:step:end这种形式的, 比如1:2:7是从1开始, 每次加2,到小于等于7为止的范围, 结果相当于[1,3,5,7],但不需要实际存储这些值。

v1[1:2:7] |> show
## [2, 5, 11, 17]

注意在倒数时,一定要指定步长-1, 像5:1这样的写法在其它编程语言中(如R)可能识别为倒数的1到5, 但是在Julia中会识别为长度为0的范围, 必须写成5:-1:1才有效地倒数。如:

v1[end:-1:begin] |> show
## [17, 13, 11, 7, 5, 3, 2]

为了获得一个向量的倒序,也可用函数reverse:

reverse(v1) |> show
## [17, 13, 11, 7, 5, 3, 2]

可以用仅有冒号作为下标, 这时表示包含所有元素的子集。 如v1[:]取出所有元素。

4.1.1.4 数组作为下标

向量的下标也可以是一个下标数组,如

v1 = [2, 3, 5, 7, 11, 13, 17]
v1[[1, 3, 5]]
    3-element Vector{Int64}:
      2
      5
     11

取出的多个元素可以修改,可以用.=运算符赋值为同一个标量,如:

v1[1:3] .= 0; 
@show v1;
## v1 = [0, 0, 0, 7, 11, 13, 17]

也可以分别赋值,如

v1[[1, 3, 5]] .= [101, 303, 505]; 
@show v1;
## v1 = [101, 0, 303, 7, 505, 13, 17]

其中.=称为广播赋值, 是将对应的元素进行赋值。

4.1.1.5 用条件作为下标

取数组子集的方括号内可以写一个等长的布尔型数组, 用来取出真值对应的项。如:

v1 = [2, 3, 5, 7, 11, 13, 17]
v1[v1 .<= 10]
4-element Vector{Int64}:
 2
 3
 5
 7

其中“.<=”表示每个元素与标量10比较。 也可以用其它向量作为条件。

数组比较的结果是BitVector:

v1 = [2, 3, 5, 7, 11, 13, 17]
show(v1 .<= 10)
## Bool[1, 1, 1, 1, 0, 0, 0]
v1 .<= 10
7-element BitVector:
 1
 1
 1
 1
 0
 0
 0

这实际是元素为布尔型(Bool)的一维数组的压缩存储版本, 每个元素仅需要一个比特存储。 可以用Array(x)将BitArray类型转换为布尔型数组, 用BitArray(x)将布尔型数组转换为压缩存储的BitArray类型。

4.1.1.6 下标运算看成函数

设长度为\(n\)的一维数组\(v\)的元素记为\(v_i, i=1,2,\ldots,n\)。 则可以将v[i]看成是如下的函数调用:

function fv(i)
    v[i]
end

这个函数的定义域是\(\{1,2,\ldots,n \}\), 值域是\(\{v_1, v_2, \dots, v_n \}\), 将\(i\)映射到\(v_i\)。 而且,可以计算v[I], 其中I是若干个下标组成的数组, 求得这几个下标映射到的\(v_i\)的值组成的向量\(\{ v_i, i \in I \}\)

如果长度为\(n\)的向量\(P\)中包含了\(1:n\)的一个排列, 则\(P\)定义了\(1:n\)的一个次序。 如:

v1 = [2, 3, 5, 7, 11, 13, 17]
P = [5, 1, 3, 4, 2, 6, 7]
@show w1 = v1[P];
## w1 = v1[P] = [11, 2, 5, 7, 3, 13, 17]

这样的\(P\)看成是一个映射, 则它是\(\{1,2,\ldots,n \}\)到自身的一个一一对应, 所以存在逆映射, 记逆映射为向量\(Q\), 则\(P[Q] = 1:n\), 如果\(w = v[P]\), 则\(w[Q] = v\)。 用sortperm(P)可以得到P的逆映射Q。如:

Q = sortperm(P)
@show P[Q];
## P[Q] = [1, 2, 3, 4, 5, 6, 7]
@show w1[Q];
## w1[Q] = [2, 3, 5, 7, 11, 13, 17]

实际上, 对长度为n的向量xsortperm(x)返回1:n的一个排列Q, 使得x[Q]从小到大排列。

4.1.2 数组类型

当一维数组元素都是整数时, 显示其类型为“Array{Int64}”, 常用的还有“Array{Float64}”, “Array{String}”, “Array{Any}”等, “Any”是Julia语言类型系统的根类型, 相应的数组可以容纳任何Julia对象作为元素。

如果数组元素都是基本类型如Float64, 则不允许给元素赋值为其它类型,如:

vf = [1.2, 2.5, 3.6]
vf[2] = "abc"
## MethodError: Cannot `convert` an object of type String to an object of type Float64
## Closest candidates are: ......

eltype()求元素类型,如:

eltype(vf)
## Float64

4.1.3 向量初始化

zeros(n)可以生成元素类型为Float64、元素值为0、长度为n的向量,如

zeros(3)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

zeros(Int64, 3)可以生成指定类型的(这里是Int64)初始化向量。如

zeros(Int64, 3)
3-element Vector{Int64}:
 0
 0
 0

Vector{Float64}(undef, n)可以生成元素类型为Float64的长度为n的向量, 元素值未初始化,如

Vector{Float64}(undef, 3)
3-element Vector{Float64}:
 1.40319995e-315
 1.40320011e-315
 1.40320027e-315

类似可以生成其它元素类型的元素值未初始化向量,如

Vector{Int}(undef, 3)
3-element Vector{Int64}:
 1800654928
          0
 1800759648

用这样的办法为向量分配存储空间后可以随后再填入元素值。

4.1.4 向量与标量的运算

向量与一个标量作四则运算, 将运算符前面加句点“.”:

    .+   .-   .*   ./   .^

表示向量的每个元素分别与该标量作四则运算, 结果仍是向量。 这种操作称为广播。 如

v2 = [2, 3, 5]
@show v2 .+ 100;
## v2 .+ 100 = [102, 103, 105]
@show 100 .- v2;
## 100 .- v2 = [98, 97, 95]
@show v2 .* 2;
## v2 .* 2 = [4, 6, 10]
@show v2 ./ 10;
## v2 ./ 10 = [0.2, 0.3, 0.5]
@show v2 .^ 2;
## v2 .^ 2 = [4, 9, 25]

为了给多个元素同时赋值为某个数,可以用“.=”操作符,如:

v2[1:2] .= 0; 
@show v2;
## v2 = [0, 0, 5]

4.1.5 向量与向量的四则运算

两个等长的向量之间作加点的四则运算,表示对应元素作相应的运算。如

v1 = [1, 3, 4, 9, 13]
v3 = [2, 5, 6, 7, 10]
@show v1 .+ v3;
## v1 .+ v3 = [3, 8, 10, 16, 23]

向量加减法也可以用不加点的运算符, 但加点版本效率略高:

@show v1 + v3;
## v1 + v3 = [3, 8, 10, 16, 23]

对应元素间乘除必须用加点的“.*”和“./”:

@show v1 .* v3;
## v1 .* v3 = [2, 15, 24, 63, 130]
@show v1 ./ v3;
## v1 ./ v3 = [0.5, 0.6, 0.6666666666666666, 
##   1.2857142857142858, 1.3]

4.1.6 向量的比较运算

两个标量之间可以进行如下的比较运算:

    ==  !==  <  <=  >  >=

结果为逻辑型(Bool),如:

5 > 3
## true

向量每个元素与标量之间、两个向量对应元素之间比较, 只要在前面的比较运算符前面增加句点:

    .==  .!=  .<  .<=  .>  .>=

运算结果是BitVector,即用1、0表示的压缩存储的真、假值。

标量比较结果之间可以用&&表示“同时成立”, ||表示“至少其中之一成立”。 布尔型标量与向量、向量之间可以用.&表示元素间“与”, .|表示元素间“或”。

例:

[1,3,5,7] .< [2,3,5,7]
4-element BitVector:
 1
 0
 0
 0

要判断向量x的每个元素是否属于向量y, 不能直接写成x .∈ y, 因为这会理解成对应元素的属于关系, 要写成x .∈ Ref(y)。如:

[1, 3] .∈ [2,3,5,7,11]
## DimensionMismatch("arrays could not be broadcast to a common size; 
## got a dimension with lengths 2 and 5")
## ......
[1, 3] .∈ Ref([2,3,5,7,11])
2-element BitVector:
 0
 1

4.1.7 向量的循环遍历

可以用for循环和eachindex()对向量的每个元素下标遍历访问。 使用in关键字,格式如

v1 = [1, 3, 4, 9, 13]
for i in eachindex(v1)
  println("v1[", i, "] = ", v1[i])
end
v1[1] = 1
v1[2] = 3
v1[3] = 4
v1[4] = 9
v1[5] = 13

这里i是循环产生的向量下标。

也可以直接写下标范围,使用“=”或者in,如:

for i = 1:length(v1)
    println("v1[", i, "] = ", v1[i])
end

因为Julia用OffsetArrays支持不从1开始的下标, 所以保险起见, 可以用firstindex(v)获取向量v的第一个下标, 用lastindex(v)获取v的最后一个下标,如:

for i = firstindex(v1):lastindex(v1)
    println("v1[", i, "] = ", v1[i])
end

也可以不利用下标而是直接对元素遍历,如

for xi in v1
    println(xi)
end
1
3
4
9
13

还可以用enumerate同时对向量下标和元素值遍历,如:

for (i, xi) in enumerate(v1)
    println("v1[", i, "] = ", xi)
end
v1[1] = 1
v1[2] = 3
v1[3] = 4
v1[4] = 9
v1[5] = 13

4.1.8 列表推导(Comprehension)

在方括号内,写[f(x) for x in collection], 可以对某个容器所有元素遍历, 并进行转换。 这称为列表推导(comprehension), 方括号内的部分称为一个生成器(generator)。 例如:

[x*x for x in [2,3,5,7]]
4-element Vector{Int64}:
  4
  9
 25
 49

还可以在末尾加上一个if cond(x)部分, 结果中仅包含满足条件的元素,如:

[x*x for x in 1:10 if x % 3 == 1] |> show
## [1, 16, 49, 100]

4.1.9 向量化运算

其它的支持向量、矩阵运算的编程语言如MATLAB、R、Python numpy都支持两个向量之间作x+y, x-y, x*y, x/y表示元素之间的运算, 以及用f(x)表示对x的每个元素调用f()变换。 Julia仅支持x+yx-y, 这些向量化运算, 有四种Julia做法:

  1. 用循环,如: julia y = similar(x) for i in eachindex(x) y[i] = f(x[i]) end

  2. 用加点广播形式,如: julia y = f.(x)

  3. 用列表推导(comprehension)形式,如: julia y = [f(xi) for xi in x]

  4. map函数,如: julia y = map(f, x)

这些做法有时可以互相替代, 有时则不方便替代。 显然加点广播形式比较易用, 但受限制也比较大,在两层循环时加点广播形式不易表达; 列表推导形式最灵活; map()往往在加点广播形式不能使用时可以替代。

4.1.10 向量的输出

在脚本文件中, 为了显示某个向量, 可以用show()函数或@show宏。 show()函数以简化格式显示向量, 不返回值(即返回nothing); @show宏已简化形式显示带有表达式提示的值, 并返回表达式结果。

比如:

v2 = [1.5, 3, 4, 9.12]
show(v2)
## [1.5, 3.0, 4.0, 9.12]
@show v2;
## v2 = [1.5, 3.0, 4.0, 9.12]

为了将向量v2按文本格式保存到文件“tmp1.txt”中, 可用:

using DelimitedFiles
writedlm("tmp1.txt", v2, ' ')

结果文件中每个数占一行。

4.1.11 向量的输入

假设文件“vecstore.txt”中包含如下的内容:

1.2 -5  3.6
7.8 9.12 4.11

可以用如下代码将文件中的数据读入到一个向量v4中:

using DelimitedFiles
v4 = readdlm("vecstore.txt")[:];

4.1.12 向量的绑定

由于Julia的变量仅仅是向实际存储空间的引用(reference), 或称绑定(binding), 所以两个变量可以引用(绑定)到同一个向量的存储空间, 修改了其中一个变量的元素值,则另一个变量的元素也被修改了。 用copy(x)可以生成x的一个浅层的副本, 如果x的元素有可变的复合类型数据, 这个浅层副本仍与原来的x有联系; 用deepcopy(x)可以生成x的一个深层副本, 一般不再与原来的x产生关联。 参见节2.1.4的说明。

4.1.13 向量的计算用函数

x是向量, sum(x)求各个元素的和, prod(x)求各个元素的乘积。

rand(n)可以用来生成n个标准均匀分布的随机数, 结果为双精度向量。 randn(n)可以用来生成n个标准正态分布的随机数, 结果为双精度向量。

maximum(v)求最大值, minimum(v)求最小值, argmax(v)求最大值首次出现的下标, argmin(v)求最大值首次出现的下标。 findmax(v)findmin(v)则返回最值和相应的首次出现的下标。

对于布尔型数组, all(v)判断所有元素为真, any(v)判断存在真值元素。

4.1.14 排序

sort(x)返回排序结果, sort(x, rev=true)返回逆序排序结果, sort!(x)原地排序。 如:

sort([11, 13, 12]) |> show
## [11, 12, 13]

sortperm(x)返回为了将x从小到大排序所应该使用的下标序列, 如:

sortperm([11, 13, 12]) |> show
## [1, 3, 2]

在有多个等长向量需要按照其中一个的次序排序时, 可以先产生这个下标序列, 然后用这个下标序列对每个要排序的向量取子集, 即可实现同时排序。

4.1.15 向量化函数

许多现代的数据分析语言, 如Python, Matlab, R等都存在循环的效率比编译代码低一两个数量级的问题, 在这些语言中, 如果将对向量和矩阵元素的操作向量化, 即以向量和矩阵整体来执行计算, 就可以利用语言内建的向量化计算获得与编译代码相近的执行效率。

Julia语言依靠其LLVM动态编译功能, 对向量和矩阵元素循环时不损失效率, 用显式循环处理向量、矩阵与向量化做法效率相近, 有时显式循环效率更高。 但是,向量化计算的程序代码更简洁, 比如上面的两个向量之间的四则运算的加点格式。

Julia中的函数, 包括自定义函数, 如果可以对单个标量执行, 将函数名加后缀句点后, 就可以变成向量化版本, 对向量和矩阵执行。 这称为广播。 如

sqrt.([1,2,3])
3-element Array{Float64,1}:
 1.0
 1.4142135623730951
 1.7320508075688772

这种向量化对于多个自变量的函数也成立。 通过编译优化, 可以达到专用的以向量作为输入的函数的效率, 如果同一个语句中有多个加点运算, 编译时能够将其尽可能地合并到同一循环中。

4.2 矩阵

前面讲的向量当元素类型相同时可以看作一维数组,不区分行向量还是列向量, 在参与矩阵运算时看作列向量。

矩阵是二维数组,有两个下标:行下标和列下标。

数组(Array)是Julia中的数据类型, 有一维、二维、多维等, 区别在于引用一个元素时所用下标个数, 数组中的元素一般都属于相同的基本类型, 比如, 元素类型都是Int64, 都是Float64, 都是String, 等等。

为了在程序中直接输入一个矩阵, 可以在方括号内两个同行的元素之间用空格分隔, 两行之间用分号分隔,如

A1 = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

注意结果显示这是一个2×3 Array{Int64,2}, 即一个2行3列的元素都是Int64类型的二维数组。

对矩阵xsize(x,1)返回行数,size(x,2)返回列数, size(x)返回行数和列数的二元组。如

@show size(A1, 1), size(A1, 2);
## (size(A1, 1), size(A1, 2)) = (2, 3)

输入行向量的例子:

[1 2 3]
1×3 Array{Int64,2}:
 1  2  3

Julia将一维向量看作列向量。所以,如下的用分号分隔的输入列向量的程序得到的是一维数组而不是二维数组:

[1; 2; 3]
3-element Array{Int64,1}:
 1
 2
 3
hcat([1,2,3])
3×1 Matrix{Int64}:
 1
 2
 3

hcat(v1, v2)将两个向量左右合并为一个矩阵。 仅输入一个向量, 就产生一个列向量。

vcat(v1, v2)将两个向量上下合并为一个矩阵, 如果仅输入一个向量, 结果仍为一维数组(看成列向量)。

对于数值型的一维数组(向量)x, 可以用x'返回x的行向量形式。

4.2.1 矩阵下标

A是矩阵,则A[i,j]表示A的第i行第j列元素,如:

A1 = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
A1[2,3]
## 6

给元素赋值可以在矩阵中修改元素值,如:

A1[2,3] = -6; A1
2×3 Matrix{Int64}:
 1  2   3
 4  5  -6

4.2.1.1 矩阵列和行

A是矩阵, 则A[:, j]表示A的第j列元素组成的向量(一维数组),如

A1 = [1 2 3; 4 5 6]
A1[:, 2]
2-element Vector{Int64}:
 2
 5

A[i, :]表示A的第i行元素组成的向量(一维数组),如

A1[2, :]
3-element Vector{Int64}:
 4
 5
 6

取出后的列或者行都是一维数组, 不分行向量和列向量。

4.2.1.2 子矩阵

如果A是矩阵,IJ是范围或者向量, 则A[I,J]表示A的行号在I中的行与列号在J中的列交叉所得的子矩阵,如

A1 = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
A1[1:2, 2:3]
2×2 Matrix{Int64}:
 2  3
 5  6
A1[[1], [1,3]]
1×2 Matrix{Int64}:
 1  3

用冒号“:”作为行下标或列下标表示取该维的全部下标。

给子矩阵赋值为一个同样大小的子矩阵给对应元素赋值,如

A1[1:2, 2:3] = [102 103; 202 203]; A1
2×3 Matrix{Int64}:
 1  102  103
 4  202  203

4.2.1.3 散列的元素

如何访问多个不构成子矩阵的元素?

CartesianIndex(i, j)可以构成一个笛卡尔下标, 访问单个元素:

B1 = [10*i + j for i=1:3, j=1:3]; display(B1)
ic = CartesianIndex(2, 3)
@show B1[ic];
3×3 Matrix{Int64}:
 11  12  13
 21  22  23
 31  32  33


B1[ic] = 23

多个元素散列情形:

ics = CartesianIndex.([1,2,2,3], [1,3,1,2])
4-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(2, 3)
 CartesianIndex(2, 1)
 CartesianIndex(3, 2)
B1[ics]
4-element Vector{Int64}:
 11
 23
 21
 32

CartesianIndices(M)返回与M形状相同的由笛卡尔下标组成的矩阵:

CartesianIndices(B1)
3×3 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
 CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)

可以将每个元素下标存在矩阵的行中:

ijmat = [1 1; 2 3; 2 1; 3 2]
4×2 Matrix{Int64}:
 1  1
 2  3
 2  1
 3  2
[B1[i, j] for (i, j) in eachrow(ijmat)]
4-element Vector{Int64}:
 11
 23
 21
 32

在R语言中可以直接写成B1[ijmat], 但Julia中不支持, B1[ijmat]的结果是将ijmat中的下标看成是B1的按列拉直后一元数组下标, 结果形状与ijmat相同。

B1[ijmat]
4×2 Matrix{Int64}:
 11  11
 21  31
 21  11
 31  21

4.2.2 矩阵初始化

zeros(m, n)可以生成元素类型为Float64、元素值为0的\(m\times n\)矩阵,如

zeros(2, 3)
2×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0

生成全是1的矩阵:

ones(2, 3)
2×3 Array{Float64,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0

zeros()ones()中可以指定元素类型,如:

ones(Int64, 2, 3)
2×3 Array{Int64,2}:
 1  1  1
 1  1  1
ones(String, 2, 3)
2×3 Array{String,2}:
 ""  ""  ""
 ""  ""  ""

fill(x, n, m)生成值都等于x的n行m列矩阵, 元素类型同于xfill(x, n)则生成一维数组。

在用方括号格式生成向量或矩阵时, 可以在方括号前面写元素类型名称, 要求生成某种类型的数组。如:

Float64[1,3,5]
3-element Array{Float64,1}:
 1.0
 3.0
 5.0

在上例中,如果不指定类型,结果将自动解析为Int64类型, 因为给出的元素值都是整数。

又如:

Float64[1 3 5; 2 4 6]
2×3 Array{Float64,2}:
 1.0  3.0  5.0
 2.0  4.0  6.0

Array{Float64}(undef, m, n)可以生成元素类型为Float64的\(m \times n\)矩阵, 元素值未初始化,如

Array{Float64}(undef, 2, 3)
2×3 Array{Float64,2}:
 1.65107e-315  6.93238e-316  6.93611e-316
 1.14941e-315  1.68873e-315  0.0

类似可以生成其它元素类型的矩阵,如

Array{Int64}(undef, 2, 3)
2×3 Array{Int64,2}:
 183993264  140312768  140313072
 140312688  140345440          0

用这样的办法先为矩阵分配存储空间,然后可以再填入元素值。

Array{T}(undef, m, n, k[, ...])这种格式可以生成元素类型为T, 初始值无定义,维数和每一维的下标长度由undef后面的参数给出, 如Array{Char}(undef, 99)是元素为字符的长度为99的字符数组, 元素初值无定义。 Array{UInt8}(undef, 3, 2, 5)是元素类型为UInt8, \(3 \times 2 \times 5\)形状的三维数组, 可以看成是5个\(3 \times 2\)矩阵, 元素初值无定义。

若x是数组, similar(x)生成元素类型和大小都与x相同, 但是元素未初始化的数组。

4.2.3 矩阵合并

A, B, C是矩阵(或向量), [A B C]将三个矩阵横向合并, 这等同于hcat(A, B, C)[A; B; C]将三个矩阵纵向合并, 这等同于vcat(A, B, C)

reshape(x, m, n)将数组x的元素填入为\(m \times n\)的矩阵, 填入时按第一列、第二列、……这样的次序,称为列次序。 如:

reshape(101:106, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 101  103  105
 102  104  106

4.2.4 矩阵元素遍历

可以对行下标和列下标分别循环, Julia矩阵元素的存储是列次序, 即第一列在最前面,第二列次之,……。 因为存储时行下标变化最快, 所以两重循环时最好令列下标在外层, 行下标在内层。如:

A1 = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
for j = 1:size(A1,2), i = 1:size(A1,1)
  println("A1[", i, ", ", j, "] = ", A1[i, j])
end
A1[1, 1] = 1
A1[2, 1] = 4
A1[1, 2] = 2
A1[2, 2] = 5
A1[1, 3] = 3
A1[2, 3] = 6

这相当于:

for j = 1:size(A1,2)
    for i = 1:size(A1,1)
      println("A1[", i, ", ", j, "] = ", A1[i, j])
    end
end

如果不考虑效率问题,也可以逐行遍历如下:

for i = 1:size(A1,1), j = 1:size(A1,2)
  println("A1[", i, ", ", j, "] = ", A1[i, j])
end
A1[1, 1] = 1
A1[1, 2] = 2
A1[1, 3] = 3
A1[2, 1] = 4
A1[2, 2] = 5
A1[2, 3] = 6

另一种办法是用类似于向量遍历的方法,如:

for i in eachindex(A1)
  println("A1[", i, "] = ", A1[i])
end
A1[1] = 1
A1[2] = 4
A1[3] = 2
A1[4] = 5
A1[5] = 3
A1[6] = 6

从上例可以看出矩阵是按列存储的。

因为Julia允许数组在OffsetArray库支持下不以1为下标起点, 所以可以用axes(A,1)获取合法行下标的迭代器, 用axes(A,2)获取列下标的迭代器, 按列次序遍历矩阵的兼容性更强的写法是:

for j = axes(A1,2), i = axes(A1,1)
  println("A1[", i, ", ", j, "] = ", A1[i, j])
end
A1[1, 1] = 1
A1[2, 1] = 4
A1[1, 2] = 2
A1[2, 2] = 5
A1[1, 3] = 3
A1[2, 3] = 6

4.2.5 矩阵读写

设当前目录中文件“vecstore.txt”中包含如下内容:

1.2 -5  3.6
7.8 9.12 4.11

将文件中的内容每一行看作矩阵的一行, 文件中保存了一个\(2 \times 3\)矩阵。 读入方法如下:

using DelimitedFiles
Ain = readdlm("vecstore.txt"); Ain

考虑上面的Ain矩阵,为了将其按文本文件格式保存到“tmp2.txt”中, 用如下程序:

writedlm("tmp2.txt", Ain, ' ')

4.2.6 矩阵与标量的四则运算

矩阵与一个标量之间用加点的四则运算符号进行运算, 与向量和标量之间的运算类似, 表示矩阵的每个元素和该变量的四则运算, 结果仍为矩阵。 这称为广播。 如

A1 = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
A1 .+ 100
2×3 Matrix{Int64}:
 101  102  103
 104  105  106
100 .- A1
2×3 Matrix{Int64}:
 99  98  97
 96  95  94
A1 .* 2
2×3 Matrix{Int64}:
 2   4   6
 8  10  12
A1 ./ 10
2×3 Matrix{Float64}:
 0.1  0.2  0.3
 0.4  0.5  0.6
A1 .^ 2
2×3 Matrix{Int64}:
  1   4   9
 16  25  36

4.2.7 两个矩阵之间的四则运算

两个同样大小的矩阵之间用加点的四则运算符号进行运算, 表示两个矩阵的对应元素的运算。如

A2 = A1 .* 100
2×3 Matrix{Int64}:
 100  200  300
 400  500  600
A1 .+ A2
2×3 Matrix{Int64}:
 101  202  303
 404  505  606
A2 .- A1
2×3 Matrix{Int64}:
  99  198  297
 396  495  594
A1 .* A2
2×3 Matrix{Int64}:
  100   400   900
 1600  2500  3600
A2 ./ A1
2×3 Matrix{Float64}:
 100.0  100.0  100.0
 100.0  100.0  100.0

矩阵加减也可以用不加点的运算符:

A1 + A2
2×3 Matrix{Int64}:
 101  202  303
 404  505  606

A1 - A2也可以。

但是, A1 * A2表示矩阵乘法, A1 / A2表示解方程组或者最小二乘, 都不是对应元素间的运算。

4.2.8 矩阵的广播运算

设矩阵\(A\)\(m \times n\), 矩阵\(B\)\(1 \times n\), 则“A .op B”可以将\(B\)的每列的同一个数重复地与\(A\)的对应列的每一个数进行运算。如:

A1 = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
B1 = [10 20 30]
1×3 Matrix{Int64}:
 10  20  30
A1 .+ B1
2×3 Array{Int64,2}:
 11  22  33
 14  25  36

类似地,若\(B\)\(m \times 1\),保存为一维数组或仅1列的矩阵, 则\(B\)的每行的同一个数可以与\(A\)的对应行的每个元素进行运算。如:

B2 = [100 200]'
2×1 adjoint(::Matrix{Int64}) with eltype Int64:
 100
 200
A1 .+ B2
2×3 Matrix{Int64}:
 101  102  103
 204  205  206

一维数组在与矩阵进行广播运算时可以看成是列向量,如:

B3 = [1000, 2000]
2-element Vector{Int64}:
 1000
 2000
A1 .+ B3
2×3 Matrix{Int64}:
 1001  1002  1003
 2004  2005  2006

4.2.9 对每行或每列的操作

矩阵经常需要对每行或每列进行同一的操作, 比如,计算行和或者列和。 设f(x)输入向量x返回标量结果, 用mapslices(f, M, dims=1)对矩阵M的第一下标用f函数进行汇总, 结果对应于每一个第二下标的值, 这等价于对M的每一列计算了结果。

mapslices(f, M, dims=2)对矩阵M的第二下标的变化方向应用f函数进行汇总, 等价于对M的每一行计算了结果, 结果是列向量, 每个元素是原来一行的汇总结果。 如:

A1=[1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
mapslices(sum, A1, dims=2)
2×1 Matrix{Int64}:
  6
 15

mapslices(f, M, dims=1)将函数f应用到第一下标变化方向上进行汇总, 结果相当于对每列计算了结果, 结果是一个行向量, 每个元素是一列的汇总结果:

mapslices(sum, A1, dims=1)
1×3 Matrix{Int64}:
 5  7  9

这种指定对哪一个下标进行汇总的方法不太直观, 但是当数组是更多维时, 能够清楚地表达汇总计算是针对哪一个下标进行汇总。

对矩阵, 也可以用eachrow()获得每行的迭代器, 用eachcol()获得每列的迭代器, 然后配合列表推导(list comprehension)进行汇总。如:

[sum(x) for x in eachrow(A1)]
2-element Vector{Int64}:
  6
 15

注意上面的结果是一维数组而不是列信号, 与mapslices(sum, A1, dims=2)的结果数据类型有所区别。

对每列计算结果:

[sum(x) for x in eachcol(A1)]
3-element Vector{Int64}:
 5
 7
 9

4.2.10 矩阵乘法

A * B表示矩阵乘法。 如

A3 = [11 12; 21 22]
2×2 Matrix{Int64}:
 11  12
 21  22
A3 * A1
2×3 Matrix{Int64}:
  59   82  105
 109  152  195

一个矩阵与一个向量(一维数组)作矩阵乘法, 向量自动变成列向量, 结果也是一维数组而不是矩阵,如:

A3 * [1, -1]
2-element Array{Int64,1}:
 -1
 -1

行向量可以直接表示成方括号内多个数值之间用空格分隔的格式,如

[1,  -1] * [1  -1]
2×2 Matrix{Int64}:
  1  -1
 -1   1

又如

[1  -1] * A3
1×2 Matrix{Int64}:
 -10  -10

注意结果是\(1 \times 2\)矩阵,即行向量,而不是向量。 向量是一维数组,行向量是二维数组。

从以上例子可以看出在矩阵运算中向量可以看成是列向量, 矩阵乘法结果如果是列向量,也会表示成向量(一维数组)。

矩阵乘法即使结果仅有一个值, 其类型也是长度为1的一维数组而不是标量, 如:

[1 2] * [-1, 1]
1-element Vector{Int64}:
 1

4.2.11 矩阵转置和旋转

A’表示矩阵A的共轭转置,对实值矩阵就是转置。 如

A1 = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
A1'
3×2 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 1  4
 2  5
 3  6

两个向量x和y的内积用LinearAlgebra包的dot(x, y)表示,结果为一个标量。如

import LinearAlgebra
LinearAlgebra.dot([1, 2], [-1, 1])
## 1

rot180()返回旋转180度的结果, rotl90()返回向左(逆时针)旋转90度的结果, rotr90()返回向右(顺时针)旋转90度的结果。

4.2.12 矩阵求逆和解线性方程组

inv(A)表示\(A^{-1}\)。如

A4 = [1 3; 3 1]
2×2 Array{Int64,2}:
 1  3
 3  1
inv(A4)
2×2 Array{Float64,2}:
 -0.125   0.375
  0.375  -0.125

A \ B表示\(A^{-1} B\), 当B是向量或者列向量时, 就是求解线性方程组\(A x = B\)中的\(x\)。 如

A4 \ [-2, 2]
2-element Array{Float64,1}:
  1.0
 -1.0

Julia提供了线性代数计算所需的一些函数, 有些在Base包中,有些则在LinearAlgebra包中。 参见LinearAlgebra包的文档。

4.2.13 其它线性代数功能

线性代数计算有关的功能通常需要用using LinearAlgebra引入线性代数包。

  • I(5): 生成5阶单位阵。
  • I: 用在矩阵乘法中,不需要指定阶数,而从两边的矩阵获得阶数。
  • size(M): 返回行、列数。
  • transpose(M): 返回仅转置的矩阵。
  • M'或者adjoint(M): 返回共轭转置,对实矩阵即转置。
  • tr(M):迹。
  • det(M): 行列式。
  • eigvalues(M): 特征值。
  • eigvecs(M): 特征向量。
  • pinv(M): 加号广义逆(Moore-Penrose广义逆)。
  • svd(M): 分解\(M = U \Sigma V^T\),返回\(U\), \(\Sigma\)的对角线元素, \(V\)

4.2.14 下标位移数组

向量、矩阵、多维数组的下标都是从1开始的。 在有些编程问题中, 从其它的下标开始会更直观。

比如, 设向量\(x\)为一个长度为\(n\)的时间序列, 要计算如下的双侧滤波(卷积): \[ y_t = \sum_{k=-1}^1 a_k x_{t+k}, \ k=2,3,\ldots,n-1 . \]\(a\)为保存了\(a_{-1}, a_0, a_1\)的数组, 如果能直接用下标\(-1,0,1\)\(y_t\)就可以写成[sum(a[k]*x[t+k] for k=-1:1) for t=2:(length(x)-1)]

OffsetArrays包提供了向量下标位移功能。如:

using OffsetArrays
let
    a = [-0.5, 1, -0.5]
    a = OffsetArray(a, -1:1)
    x = [2,3,5,7,11,13]
    [sum(a[k]*x[t+k] for k=-1:1) for t=2:(length(x)-1)]
end
4-element Vector{Float64}:
 -0.5
  0.0
 -1.0
  1.0

对矩阵也可以如此定义下标位移, 如:

A = [10*i + j for i=1:3, j=1:3]
3×3 Matrix{Int64}:
 11  12  13
 21  22  23
 31  32  33
B = OffsetArray(A, -1:1, 0:2)
B[0,0]
## 21

4.2.15 特殊矩阵

求矩阵的对角阵:

using LinearAlgebra
A = [1 2 3; 4 5 6; 7 8 9]
Diagonal(A)
3×3 Diagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  5  ⋅
 ⋅  ⋅  9

给定对角元素生成对角阵:

Diagonal([1, 5, 9])
3×3 Diagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  5  ⋅
 ⋅  ⋅  9

生成稀疏矩阵:

using LinearAlgebra, SparseArrays
A = [1 0 4; 0 0 2; 3 0 0]
sparse(A)
3×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  ⋅  4
 ⋅  ⋅  2
 3  ⋅  ⋅

稀疏矩阵也可以用sparse(I, J, V)生成, 三个输入参数的对应元素分别为行号、列号、值。如:

sparse([1,1,2,3], [1,3,3,1], [1,4,2,3])
3×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  ⋅  4
 ⋅  ⋅  2
 3  ⋅  ⋅

系数矩阵的存储格式为压缩系数列格式(compressed sparse column (CSC))。 矩阵按列保存, 保存了行数、列数, 保存了每列首个元素在存储中的序号colptr, 如果某列都是零则该列对应的序号值等于前一列的序号值; 保存了保存的每个值所在行下标(列下标从colptr中获得), 按列保存了所有非零值。 参见:

4.3 数组

4.3.1 初始化

向量(Vector{T})是一维数组, 矩阵(Matrix{T})是二维数组, 未初始化的d维数组可以用Array{T,d}(undef, dims)生成。

Array{Int,3}(undef, (2,3,2))
2×3×2 Array{Int64, 3}:
[:, :, 1] =
         0  286486736  282133312
 282133312          0  281952336

[:, :, 2] =
         0  281952336  281952336
 282133312          0  281952336

初始化的三维数组:

x = reshape(Vector{Int}(1:12), (2,3,2))
2×3×2 Array{Int64, 3}:
[:, :, 1] =
 1  3  5
 2  4  6

[:, :, 2] =
 7   9  11
 8  10  12

三维数组显示时, 将x[i,j,k]中每一个固定k所对应的矩阵分别显示, 可以将i看成行,j看成列,k看成“层”。

可以用zeros(dims), ones(dims)生成Float64类型的数组, 用fill(xi, dims)生成元素都等于xi的数组。 可以用similar(x)生成与x同类型、同大小的未初始化数组, 用similar(x, type)生成大小相同、指定类型的数组, 用similar(x, dims)生成类型相同、指定大小的数组, 用similar(x, type, dims)生成类型、大小指定的数组。如:

similar([1 2; 3 4])
2×2 Matrix{Int64}:
 247703280  284851696
 284851696  183369904
similar([1 2; 3 4], Float64)
2×2 Matrix{Float64}:
 9.05996e-316  9.05976e-316
 9.05976e-316  9.05976e-316
similar([1 2; 3 4], (4,))
4-element Vector{Int64}:
 249499536
 249499568
 342752640
 249500976
similar([1 2; 3 4], Float64, (4,))
4-element Vector{Float64}:
 1.4017172e-315
 1.4017172e-315
 1.14117172e-315
 9.0596687e-316

4.3.2 数组的基本函数

ndims(x)返回维数, 向量为一维, 矩阵为二维。

size(x)返回各维长度组成的元组, 对矩阵即行数和列数。 可以用size(x, dim)求某一维的元素个数, 如矩阵的size(x, 2)为第二下标个数即列数。 length(x)返回总的元素个数。

axes(x)返回一个对应于各维的元组, 每个元素是该维的可用的下标。 axes(x, dim)则返回指定的哪一维的可用的下标。 如:

collect.(axes(fill(0, (2,3,4))))
## ([1, 2], [1, 2, 3], [1, 2, 3, 4])

keys(x)返回x的所有下标组合构成的数组,如:

keys(zeros(2,3))
2×3 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)

eachindex(x)返回高效遍历x的下标迭代器, 一般是按一维线性列优先方式遍历的, 某些特殊数组如稀疏矩阵可能会返回多维下标形式的迭代器。 用IndexStyle(x)查询x基本的下标访问方式, 为IndexLinear()IndexCartesian()

4.3.3 数组的广播

f.(arr)对每个元素计算f值; f.(arr...)对多个数组的对应元素计算f值, 允许其中有标量,标量被重复利用。

4.3.4 合并

cat(arrs...; dims)沿dims指定的维度进行合并, 而其它没有指定的维度必须是完全同形状的。 如:

A = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
B = [11 12; 21 22]
2×2 Matrix{Int64}:
 11  12
 21  22
cat(A, B, dims = 2)
2×5 Matrix{Int64}:
 1  2  3  11  12
 4  5  6  21  22

在上例中,AB的第一维度(行)相同, 要求在第二维度(列)方向进行合并, 结果是合并了列, 结果的行维不变。 其中的dims = 2也可以写成更一般的dims = (2,)。 如果dims中有多个维度, 结果是分块对角矩阵或者高维的分块对角矩阵类似结构。

hcat对矩阵表示左右合并, 对高维数组表示沿第二维方向合并; vcat对矩阵表示上下合并, 对高维数组表示沿第一维方向合并。

4.3.5 数组函数

cumsum(A; dims)沿指定的维度计算累计的和。 对一般的迭代器也有cumsum(itr)cumprod(A; dims)沿指定的维度计算累计的乘积。 对一般的迭代器也有cumprod(itr)

accumulate(op, A; dims, [init])cumsum的推广版本, 沿指定的维度进行op要求的累计计算。

diff(A)计算差分, diff(A; dims)沿指定的维度计算差分。

repeat(x, n)将数组x重复n次, 而repeat(x, n...)将数组x的每个元素按照多个n...指定的维度重复。如:

repeat([11,12], 3)
6-element Vector{Int64}:
 11
 12
 11
 12
 11
 12
repeat([11;12], 2, 3)
4×3 Matrix{Int64}:
 11  11  11
 12  12  12
 11  11  11
 12  12  12

在上例中,先将列向量[11;12]沿第一维度(行下标方向)进行了2次重复,得到[11;12;11;12]; 再将结果延第二维度(列下标方向)进行了3次重复,得到如上结果。

对矩阵,eachrow(A)给出各个行向量的迭代器, eachcol(A)给出各个列向量的迭代器。 对数组, eachslice(A; dim)给出关于dim指定的维度的迭代器, 效果是对dim指定的维度的每个下标产生一个子数组。 如:

A = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
sum.(eachrow(A))
2-element Vector{Int64}:
  6
 15
sum.(eachslice(A, dims = 1))
2-element Vector{Int64}:
  6
 15
A = reshape(1:12, (2,3,2))
2×3×2 reshape(::UnitRange{Int64}, 2, 3, 2) with eltype Int64:
[:, :, 1] =
 1  3  5
 2  4  6

[:, :, 2] =
 7   9  11
 8  10  12
first(eachslice(A, dims=3))
2×3 view(reshape(::UnitRange{Int64}, 2, 3, 2), :, :, 1) with eltype Int64:
 1  3  5
 2  4  6
sum.(eachslice(A, dims=3))
2-element Vector{Int64}:
 21
 57

4.3.6 排列组合

xn个元素, isperm(x)判断x是否为1:n的一个排列。

v1:n的一个排列, 将v看成是从\(1:n\)\(v_1 \sim v_n\)的映射, invperm(v)返回v的逆映射, 即若A = B[v], 则B = A[invperm[v]]。 如:

v = [2, 3, 1] # 1 => 2, 2 => 3, 3 => 1
3-element Vector{Int64}:
 2
 3
 1
invperm(v) # 1 => 2, 2 => 1, 3 => 2
3-element Vector{Int64}:
 3
 1
 2

p1:n的一个排列, permute!(v, p)p的次序重排v。 如:

v = [11,12,13]; p = [2,3,1]
permute!(v, p)
3-element Vector{Int64}:
 12
 13
 11
v = [11,12,13]
v[p]
3-element Vector{Int64}:
 12
 13
 11

invpermute!(v, p)则按照p的逆变换将v修改次序。

reverse(v)返回颠倒次序的向量, reverse(A; dims=d)返回沿指定维度颠倒次序的数组:

A = [1 2 3; 4 5 6]
reverse(A, dims=2)
2×3 Matrix{Int64}:
 3  2  1
 6  5  4

可以沿多个维度颠倒次序,如:

reverse(A, dims=(1,2))
2×3 Matrix{Int64}:
 6  5  4
 3  2  1

reverse(A, :)表示沿所有维度颠倒次序。

reverse!是原地颠倒的版本。

4.4 静态数组

StaticArrays包提供了一种从类型就固定了大小的数组, 比基本库中的Array类型减少了一些通用性, 但对小规模数组的效率大大提升, 常常有一个数量级的效率提示。 称为静态数组, 是不可修改类型; 也有可修改的变种。

参见:

可以用SASA_F64和初值生成静态数组,如:

using StaticArrays
x1 = SA[1,2,3]
3-element SVector{3, Int64} with indices SOneTo(3):
 1
 2
 3

这产生了类型SVector{3, Int64}, 这是静态一维数组, 特点是在类型中规定了数组大小,不能修改大小, 也不能修改存储的值。

类似地,可以生成指定类型的静态向量:

x2 = SA_F64[1,2,3] # 元素类型为Float64
3-element SVector{3, Float64} with indices SOneTo(3):
 1.0
 2.0
 3.0

静态矩阵:

m1 = SA[1 2; 3 4]
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1  2
 3  4

m1的类型是SMatrix{2, 2, Int64, 4}, 即元素类型为Int64,大小为\(2 \times 2\),有4个元素的静态矩阵。

Float64类型的静态矩阵:

m2 = SA_F64[1 2; 3 4]
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 1.0  2.0
 3.0  4.0

可以用SVector, SMatrix函数生成指定大小、指定元素类型的初始化的静态数组。 MVectorMMatrix是指定大小、指定元素类型但允许修改元素值的数组。 可以用@SVector, @SMatrix, @SArray将已有数组转换为静态数组。

4.5 生成器

对容器的循环,经常可以表达成一种称为生成器的语法, 写法为“f(x) for x in collection if cond(x)”。 如果写在[]内, 就生成数组, 写在Tuple()内, 就生成元组; 写成Dict(x => f(x) for x in collection if cond(x)), 就生成字典。 也可以用来为接受一个迭代器作为输入的函数提供输入。

4.5.1 生成向量

例如,为了生成1, 2, 3的立方的向量,可以写成:

xcube = [i^3 for i=1:3]
3-element Array{Int64,1}:
  1
  8
 27

对前5个素数作立方,可以写成

x = [2, 3, 5, 7, 11]
y = [z^3 for z in x]
@show y;
## y = [8, 27, 125, 343, 1331]

当然,这样的简单运算也可以用加点四则运算表示:

x = [2, 3, 5, 7, 11]
y = x .^ 3
@show y;
## y = [8, 27, 125, 343, 1331]

sum()函数可以接受一个遍历器, 所以可以用生成器作为输入:

sum(x^2 for x in 1:10 if x % 3 == 1)
## 166

4.5.2 用生成器生成矩阵

可以在方括号内写双层循环生成矩阵,如:

Ac = [100*i + j for i=1:2, j=1:3]
2×3 Matrix{Int64}:
 101  102  103
 201  202  203

这里i对应于行,j对应于列, 应该理解为: \[\begin{aligned} & \text{for } i \text{ in } 1:2 \\ &\qquad \text{for } j \text{ in } 1:3 \\ &\qquad\qquad A[i,j] = 100i + j \\ &\qquad \text{end} \\ & \text{end} \end{aligned}\] i是外层循环,j是内层循环。

Python也有类似做法但是次序相反。

这种格式也允许在方括号前写出要求的元素类型,如:

Float64[i*100 + j for i=1:2, j=1:3]
2×3 Matrix{Float64}:
 101.0  102.0  103.0
 201.0  202.0  203.0

如下程序返回矩阵对角化的结果:

[(i==j ? Ac[i,i] : 0) for i=1:2, j=1:3]
2×3 Matrix{Int64}:
 101    0  0
   0  202  0

4.5.3 用生成器生成字典

Dict{Int, Int}(i => i*i for i = 1:5)
Dict{Int64, Int64} with 5 entries:
  5 => 25
  4 => 16
  2 => 4
  3 => 9
  1 => 1

4.6 收集容器类的数据结构

数组、元组、字典、集合等都属于收集容器(collections), 它们支持许多共同的接口(函数), 如长度、遍历等。

4.6.1 迭代器

迭代器(Iterators)是一种接口, 许多收集容器支持这种接口, 而许多函数,如sumreduce, mapreduce, 只要输入支持迭代器接口的数据类型就可以给出正确结果, 即使这个数据类型是用户自定义的实现了迭代器接口的类型也可以。

支持迭代器的类型主要是实现了iterate函数。 iterate(iter)或者返回第一个元素和当前状态的元组, 或者返回nothing, 而iterate(iter, state)则从已有的状态出发访问下一状态, 或者返回下一个元素值以及当前状态, 或者返回nothing

比如,常用的for循环:

for i in iter
    # ......
end

实际上是用如下的iterate函数调用实现:

next = iterate(iter)
while next !== nothing
    (i, state) = next
    # ......
    next = iterate(iter, state)
end

其中的状态state与具体的实现有关, i是每一次迭代获得的值。

4.6.2 范围类型

AbstractRange是范围类型的父类型, 支持迭代器, 实际上本身就只是一个迭代器而不需要存储其中的元素。

OrdinalRange{T,S}支持等差的序列, 且类型T必须为整型等有最小间隔的类型, S代表公差类型。

StepRange{T,S}支持一般的等差序列, T为元素值的类型,S为公差类型。

UnitRange{T}是公差等于1的等差数列类型, 如1:5这样的写法属于UnitRange{Int64}

range()函数返回AbstractRange的等差数列。 range(start, stop, length)指定长度, range(start, stop; length, step)可选lengthsteprange(start; length, stop, step)输入起始值和length, stop, step中的两个; range(;start, length, stop, step)输入start, length, stop, step中的三个, 也可以仅输入start, length, stop中的两个而默认step=1range会以有理数而非浮点数方式计算数列的值。 如:

range(0, 1, 6)
## 0.0:0.2:1.0
range(0, 1, step=0.2)
## 0.0:0.2:1.0

4.6.3 一般收集容器的接口

isempty(collection)判断是否为空。 empty!(collection)清空内容。

length(collection)返回长度。

4.6.4 支持迭代器的收集容器接口

in, 判断属于关系,一般对每个元素用==来判断。 有missing时结果可能为missing

Set类型用isequal判断, 字典类型仅判断键(key)的属于关系。

不属于用

判断多个元素每一个是否属于某个收集容器, 不能简单地使用A .∈ S的形式, 因为这时广播作用针对AS都存在, 应改写为A .∈ Ref(S).∉类似,写成A .∉ Ref(S)。 如:

1  [1,2,3]
## true
0  [1,2,3]
## false
[0,1] .∈ Ref([1,2,3])
## 2-element BitVector:
##  0
##  1

eltype(collection)返回元素类型。如:

eltype(1:5)
## Int64

indexin(a, b)a中每个元素返回在b中首次匹配的位置下标, 不匹配的返回nothing,如:

indexin([1,3,5], [1,2,3])
3-element Vector{Union{Nothing, Int64}}:
 1
 3
  nothing

collect(iter)将一个迭代器转换为一个数组。 在列表推导结构中利用了这个函数将生成器结果转换为数组。 如:

collect(1:2:7)
4-element Vector{Int64}:
 1
 3
 5
 7

collect(element_type, collection)可以按指定元素类型返回转换后的数组,如:

collect(Float64, 1:2:7)
4-element Vector{Float64}:
 1.0
 3.0
 5.0
 7.0

unique(iter)返回容器中所有不同的元素。 unique(f, iter)返回容器中用f映射为不同值的元素。 unique!(x)修改x使其仅保留不同元素, unique!(f, x)修改x使其仅保留用f映射为不同值的元素。 可以用allunique(iter)判断是否元素彼此不同。

除了sum(iter; [init]), 还可以用sum(f, iter; [init])求函数值的和, 关键字参数init用来给出当iter为空时的函数值和累加初值。 prod是类似的函数,计算乘积。 cumsum计算并输出累加序列, cumprod计算并输出累乘序列。

maximum(iter)求最大值, maximum(f, itr; [init])求用f变换后的最大函数值, init提供了itr为空时的函数值。

argmax(iter)求最大值所在的下标, argmax(f, iter)f的最大值点。 findmax(f, inter)返回最大值和最大值点下标位置的元组。

最小值有类似的函数minimum, argmin, findmin等。

any(iter)判断是否至少一个元素为真, any(p, iter)用示性函数p判断至少一个元素满足条件。 all则判断是否所有元素为真或满足条件。 count(iter)计算true的个数, count(p, iter)计算满足示性函数p给出的条件的元素个数。

first(iter)求得第一个元素, first(iter, n)返回前n个元素, last(iter)求得最后一个元素, last(iter, n)返回最后n个元素。

front(tup)返回元组除去最后一个元素后构成的元组, tail(tup)返回元组除去第一个元素后构成的元组。

4.6.5 关于迭代器的方便函数

Stateful(iter)用来生成可修改的、保存状态的迭代器。

zip(iters...)可以将多个迭代器按照对应元素合并, 结果为合并元素得到的元组的迭代器。

enumerate(iter)返回关于(i, xi)的迭代器, i是序号, xi是元素值, 用如:

for (i, xi) in x:
    println("x[", i, "] = ", xi)
end

rest(iter, state)返回从状态state开始的迭代器,如:

itr = Iterators.rest([11, 13, 17, 19], 3)
collect(itr)
    2-element Vector{Int64}:
     17
     19

Iterators.countfrom(start, step)生成一个仅有开始没有结束的迭代器,如:

for x in Iterators.countfrom(11, 2)
    println(x)
    if(x >= 15) 
        break
    end
end
## 11
## 13
## 15

Iterators.take(iter, n)取出迭代器的前n个元素(不足n个时返回所有元素)的迭代器, 如:

Iterators.take((x^2 for x = 1:10), 3) |> collect
3-element Vector{Int64}:
 1
 4
 9

Iterators.takewhile(p, itr)输出满足示性函数p条件下的迭代器, 类似于while循环。 如:

Iterators.takewhile(xi -> xi <= 17, [11, 13, 17, 19]) |> collect
3-element Vector{Int64}:
 11
 13
 17

Iterators.drop(itr, n)返回去掉前n个后的迭代器。如:

Iterators.drop([11, 13, 17, 19], 2) |> collect
2-element Vector{Int64}:
 17
 19

Iterators.dropwhile(p, itr)舍弃满足示性函数p条件的前面若干个值后返回剩余值的迭代器, 如:

Iterators.dropwhile(xi -> xi <= 13, [11, 13, 17, 19]) |> collect
2-element Vector{Int64}:
 17
 19

Iterators.cycle(itr)返回对输入的迭代器无穷地重复的迭代器。

Iterators.repeated(xi)返回无穷地重复xi值的迭代器。 Iterators.repeated(x, n)返回重复xin次的迭代器。

Iterators.product(itrs...)对多个迭代器进行完全组合, 返回结果的迭代器,为多维数组形式, 每个元素为组合的元组。如:

Iterators.product([2,3,5,7], [11, 13]) |> collect
4×2 Matrix{Tuple{Int64, Int64}}:
 (2, 11)  (2, 13)
 (3, 11)  (3, 13)
 (5, 11)  (5, 13)
 (7, 11)  (7, 13)

注意这与zip不同,zip是将对应元素构成元组,如:

zip([2,3,5,7], [11, 13]) |> collect
2-element Vector{Tuple{Int64, Int64}}:
 (2, 11)
 (3, 13)

Iterators.flatten(itrs)可以将元素为迭代器的迭代器展开为单个的迭代器,如:

Iterators.flatten( ((1,2), (3, 4, 5)) ) |> collect
5-element Vector{Int64}:
 1
 2
 3
 4
 5

Iterators.partition(itr, n)将迭代器按n个一段分成若干段, 输出为迭代器。如:

Iterators.partition(1:10, 3) |> collect
4-element Vector{UnitRange{Int64}}:
 1:3
 4:6
 7:9
 10:10

Iterators.map(f, itrs...)产生一个用f作用后的迭代器, 但仅使用时才真的调用f进行计算。 Base.map(f, itrs...)则需要对元素调用f进行计算。

Iterators.filter(p, itr)产生一个用p过滤后的迭代器, 仅使用时才调用p过滤。

Iterators.reverse(iter)生成一个反向的迭代器, 与Base.reverse(iter)相比的区别是Iterators.reverse(iter)仅在使用时才计算。

Iterators.only(iter)返回iter的仅有的一个元素, 否则出错。

Iterators.peel(iter)返回第一个元素和剩余部分的迭代器。

4.6.5.1 Iterators.product例子:完全试验方案

在多因素试验设计中, 可以用Iterators.product生成完全试验的试验方案。 比如,因素A有水平A1, A2, 因素B有水平B1, B2, 因素C有水平C1, C2, 产生完全搭配的数组,如:

facs = [ ["A1", "A2"], ["B1", "B2"], ["C1", "C2"] ]
collect(Iterators.product(facs...))
2×2×2 Array{Tuple{String, String, String}, 3}:
[:, :, 1] =
 ("A1", "B1", "C1")  ("A1", "B2", "C1")
 ("A2", "B1", "C1")  ("A2", "B2", "C1")

[:, :, 2] =
 ("A1", "B1", "C2")  ("A1", "B2", "C2")
 ("A2", "B1", "C2")  ("A2", "B2", "C2")

其中facs...是将一个列表的元素作为某个函数的各个自变量的写法。 为了将这些试验组合转换为每行一个试验组合的格式, 用[:]的做法将多维数组转换为按列拉直的一维数组, 每个元素是一个试验方案组合:

collect(Iterators.product(facs...))[:]
8-element Vector{Tuple{String, String, String}}:
 ("A1", "B1", "C1")
 ("A2", "B1", "C1")
 ("A1", "B2", "C1")
 ("A2", "B2", "C1")
 ("A1", "B1", "C2")
 ("A2", "B1", "C2")
 ("A1", "B2", "C2")
 ("A2", "B2", "C2")

可以用collect.将每个试验方案从元组转换为一维数组, 得到一个一维数组的一维数组:

collect.(collect(Iterators.product(facs...))[:])
8-element Vector{Vector{String}}:
 ["A1", "B1", "C1"]
 ["A2", "B1", "C1"]
 ["A1", "B2", "C1"]
 ["A2", "B2", "C1"]
 ["A1", "B1", "C2"]
 ["A2", "B1", "C2"]
 ["A1", "B2", "C2"]
 ["A2", "B2", "C2"]

reduce, reshapevcat将上述的一维数组的一维数组转换成矩阵, 矩阵每行为一个试验方案:

Iterators.product(facs...) |>
    collect |>
    xs -> xs[:] |>
    xs -> collect.(xs) |>
    xs -> [reshape(x, (1,3)) for x in xs] |>
    xs -> reduce(vcat, xs)
8×3 Matrix{String}:
 "A1"  "B1"  "C1"
 "A2"  "B1"  "C1"
 "A1"  "B2"  "C1"
 "A2"  "B2"  "C1"
 "A1"  "B1"  "C2"
 "A2"  "B1"  "C2"
 "A1"  "B2"  "C2"
 "A2"  "B2"  "C2"

写成一个函数,用来生成多因素试验的试验方案组合, 输出一个矩阵, 每行为一个组合:

function allcomb(factors)
    numfactors = length(factors)
    Iterators.product(factors...) |>
        collect |>
        xs -> xs[:] |>
        xs -> collect.(xs) |>
        xs -> [reshape(x, (1,numfactors)) for x in xs] |>
        xs -> reduce(vcat, xs)
end 
allcomb(facs)
8×3 Matrix{String}:
 "A1"  "B1"  "C1"
 "A2"  "B1"  "C1"
 "A1"  "B2"  "C1"
 "A2"  "B2"  "C1"
 "A1"  "B1"  "C2"
 "A2"  "B1"  "C2"
 "A1"  "B2"  "C2"
 "A2"  "B2"  "C2"

4.6.6 支持下标运算的收集容器接口

A[i,j]这样的下标运算, 实际是调用了getindex(A, i, j)。 只要定义了getindex函数, 就可以支持下标运算。 一般格式为getindex(collection, key...)。 对数组, 下标是一个或多个整数; 对字典, 下标是不可修改类型, 常用字符串。

还要定义setindex!(collection, value, key...)用来实现A[i,j]=b这样的赋值操作。

firstindex(x)用来定义x[begin]firstindex(x, d)用来定义第d维下标的第一个值。 每一个下标的第一个值一般等于1, 但可以用下标位移方法定义不从1开始的下标值。

lastindex(x)用来定义x[end]lastindex(x,d)用来定义第d维下标的结束值。

4.6.7 AbstractArray的接口

必须有getindexsetindex, firstindex, lastindex

size(A)返回包含各维长度的元组。 length(A)返回总元素个数。

可支持iterate,按合并各维、列优先次序遍历。

可支持similar(A),生成与A维数相同、元素类型相同的未初始化数组。 可支持similar(A, ::Type{S}, dims::Dims)其中元素类型和维数都可以变化, 也可以有similar(A, ::Type{S}),仅改变元素类型, 或similar(A, dims::Dims),仅改变维数。

AbstractArray为父类的类型都支持按一个下标遍历或访问, 按多个下标访问。

AbstractArray是参数化类型, 继承时用MyType <: AbstractArray{T, N}的格式, T是元素类型,如Float64, N是维数,如1。

在自定义继承了AbstractArray的类型时, 类型按下标基本运算分为IndexLinear()IndexCartesian()两种, 只需要定义其中效率最高的一种, 另一种可以自动转换。 稀疏矩阵就是IndexCartesian()类型的。

4.6.8 字典的接口

Julia在标准库中提供了Dict类型, 用杂凑表(hash table)保存数据, 用isequal匹配键。

可以继承AbstractDict{S,T}定义新的字典类型。 必须实现haskey(x, key), 表示key键是否存在。 keys(x)返回所有键的迭代器,次序不保证。 values(x)返回所有值的迭代器,次序不保证,但是与keys(x)次序一致。 pairs(x)返回键和值的元组的迭代器,次序不保证。 get(x, key, default)取出key指向的值,key不存在时返回defaultget!(x, key, default),这个版本在key不存在时将key => default插入到字典中并返回该值。 get!(f, x, key)key不存在时将key => f(key)插入到字典中并返回该值。

getkey(x, key, default)返回匹配的key本身而不是指向的值, 无匹配时返回default

keytype(x)返回键的数据类型。 valtype(x)返回值的数据类型。

delete!(x, key)删除指定的键值对, 返回修改后的字典。 pop!(x, key[, default])删除key指定的键值对并返回该值, 不匹配时返回default

merge(dict1, dict2)合并两个字典。 存在相同键时以第二个字典的值为准。 merge!(dict0, dict2)用第二个字典增补修改第一个字典。

4.6.9 集合的接口

Julia标准库实现了Set, 也可以继承AbstractSet{T}实现自己的集合类型, T是元素类型。

Set([itr])从一个遍历器生成集合。

union(A, B)A ∪ B可计算两个收集容器的并集, union!(A, B)则直接修改A

intersect(A, B)A ∩ B计算交集。 intersect!(A, B)则修改A

setdiff(A, B)计算差集A\Bsetdiff!(A,B)则修改Asymdiff(A,B)计算双向差集的并集,称为对称差集。 symdiff!(A,B)则修改A

issubset(A, B)A ⊆ BB ⊇ A判断子集关系是否成立。 否定这种关系用A ⊈ BB ⊉ B判断。 输入方法:

  • : \subseteq<TAB>
  • : \supseteq<TAB>
  • : \nsubseteq<TAB>
  • : \nsupseteq<TAB>

issetequal(A, B)判断作为集合是否相等, 不计次序和重复。

isdisjoint(A,B)判断是否不交(交集为空)。

4.7 map-reduce等函数式编程

4.7.1 map

map(f, c)对迭代器c的每个元素用f变换并返回结果, 或对多个迭代器的对应元素用f变换并返回结果, 当其中某一个迭代器元素用完时遍历结束。

单个迭代器的map往往可以用生成器或者广播函数代替。 例:

map(x -> x^2, [1,3,5])
3-element Vector{Int64}:
  1
  9
 25
[x^2 for x in [1,3,5]]
3-element Vector{Int64}:
  1
  9
 25
(x -> x^2).([1,3,5])
3-element Vector{Int64}:
  1
  9
 25

有多个迭代器时map不要求长度相等, 这一点与广播不同。如:

map(+, [10,20], [1,3,5])
2-element Vector{Int64}:
 11
 23

map!(function, destination, collection...)可以将结果保存在输入的destination中。

对字典dict, map!(f, values(dict))可以用f修改dict的值。

foreach(f, c...)map的变种, 差别在于不返回结果。

map(f, collection)也可以写成

map(collection) do x
    # f(x)的函数体语句
end

如:

map(1:5) do x
    x^2
end
5-element Vector{Int64}:
  1
  4
  9
 16
 25

有时为了重复某个语句n次而不想写成for循环, 也可以写成map(1:n) do _ ... end的形式, 结果为长度\(n\)的一维向量, 如:

map(1:5) do _
    rand(11:15)
end
5-element Vector{Int64}:
 14
 12
 11
 14
 13

4.7.2 reduce

reduce(op, itr; [init])将二元的运算op推广到多元运算, 实际就是按照依次地用op进行汇总。 可以用init指定一个初始值。如:

reduce(*, 1:5) == prod(1:5)
## true
reduce(*, 1:5, init=-1)
## -120
reduce(, [[1,2], [1,3], [2,4]]) |> show
## [1, 2, 3, 4]

reduce实际上是不规定运算次序是从左到右还是从右到左的, 可以改用foldl要求从左到右计算, 用foldr要求从右到左计算。

许多常用的reduce操作以及提供了单独的函数, 如sum, prod, minimum, maximum, all, any

在实际中经常进行迭代计算, 若每次迭代的结果是一个长度为\(m\)的一维数组, 则在已知共\(n\)次迭代时, 应预先分配一个\(m \times n\)矩阵\(M\), 将第\(j\)次迭代的结果保存到\(M\)的第\(j\)列中。 但是, 迭代也经常不能预知需要进行的迭代次数, 这时, 可以将结果保存成一维数组的一维数组, 如:

result = []
while loop_cond
    new_item = step()
    push!(result, new_item)
end

其中结果new_item是长度为\(n\)的一维数组, result是一维数组的一维数组。 在结束迭代后, 应将结果保存为矩阵形式以利于进一步处理。 这时可利用reduce:

result_matrix = reduce(hcat, result)

这在各次迭代的结果之间进行横向合并, 得到一个\(n \times m\)矩阵。如:

res = [[1,2,3], [4, 5, 6], [7, 8, 9]]
reduce(hcat, res)
3×3 Matrix{Int64}:
 1  4  7
 2  5  8
 3  6  9

4.7.3 accumulate

sum将加法从二元运算扩充为多元运算, reduce将一般的满足交换律、结合律的二元运算扩充为多元运算; cumsum将加法扩充成累加运算, accumulate则加一般的满足交换律、结合律的二元运算扩充累计的运算。 accumulate(f, itrs)对运算f进行累计运算, accumulate(f, itrs, init)则增加一个计算初值。 如:

accumulate(, [ [2,11], [2,3,5], [7,11] ])
3-element Vector{Vector{Int64}}:
 [2, 11]
 [2, 11, 3, 5]
 [2, 11, 3, 5, 7]

4.7.4 mapreduce

mapreduce(f, op, itrs...; [init])等价于 reduce(op, map(f, itrs...); [init]), 但不需要保存中间结果, 效率会略高一些。 如:

mapreduce(x -> x^2, +, 1:5)
## 55

上式计算了\(\sum_{k=1}^5 k^2\)

reduce类似, mapreduce不保证汇总计算时的次序, 而mapfoldl总是从左向右计算, mapfoldr总是从右向左计算。

4.7.5 filter

filter(f, iter)用示性函数f挑选iter的元素, 返回挑选的结果。如:

filter(iseven, 1:10)
5-element Vector{Int64}:
  2
  4
  6
  8
 10

如果iter是字典, f将应用到key => value对上, 挑选满足条件的条目, 返回保留的内容构成的字典。

如果iterskipmissing输出的迭代器, 则输出非缺失、f值为真的哪些元素。

filter!(f, x)则修改x使其仅保留示性函数f值为真的元素。

4.7.6 replace

replace(A, old_new::Pair...)可以从迭代器A中用给出的一个或多个替换规则进行替换。 如:

replace([1, 3, 5, 3, 5, 7], 3 => 30) |> show
## [1, 30, 5, 30, 5, 7]
replace([1, 3, 5, 3, 5, 7], 3 => 30, 5 => -5) |> show
## [1, 30, -5, 30, -5, 7]

可以加选项count=n指定最多替换n次。

replace(f, iter; count=n)可以用f变换替换iter的元素值, 可以用关键字参数count指定最多替换多少次。 替换次数是按变换后的值不等于原始值的次数计算的。 如果都替换, 可以用map或广播。

replace!(A, old_new::Pair...; count=n)replace!(f, iter; count=n)进行原地替换, 这是map不能取代的。

4.7.7 findfirst

findfirst(p, x)可以找到满足示性函数p的条件的下标。如:

findfirst(xi -> xi==12, 11:15)
## 2

类似函数还有findlast, findall等。

对于最大值点, 可以用argmax(x), argmax(f, x)等求得最大值点下标。