5
4

More than 5 years have passed since last update.

[Julia 0.5] 配列添字を0から数える

Last updated at Posted at 2016-12-23

Julia Advent Calendar 2016 です。12月2日の埋め草として 12月23日に投稿しました。一人で3つも書くのは反則かしらん。

 Julia では、配列の添字を 1 から数えますが (1-based array index)、Julia 0.5 では、1 以外から添字を数える機能が試験的に実装されました。該当のドキュメント Arrays with custom indices には、この機能の利用上の注意しか記載されていません。何らかのパッケージを通じて使うようです。
 ここでは、OffsetArray パッケージを紹介します。Qiita記事 ラプラス方程式 内の Fortran 2008 プログラムを Julia に移植してみました。

下の Jupyter notebook: https://gist.github.com/anonymous/7437f799d132fb15e4cf0e086d5c16a9

using PyPlot
using OffsetArrays
n = 100

# real :: mesh(0:n + 1, 0:n + 1)
# mesh(0:n + 1, 0:n + 1) # Good
mesh0 = Array( Float64, n+2, n+2) 
mesh = OffsetArray(mesh0, 0:n+1, 0:n+1)

mask = falses( (n,n) )

## define shape : eccentric pipe   
#  forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mask(i, j) = .true.
for i = 1:n, j = 1:n
    if (i - 51)^2 + (j - 51)^2 < 50^2
        mask[i, j] = true
    end
end

# forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 71)**2 < 10**2) mask(i, j) = .false.
for i = 1:n, j = 1:n
    if (i - 51)^2 + (j - 71)^2 < 10^2
        mask[i, j] = false
    end
end

# mesh=100.0 
fill!(mesh, 100.0)

# forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mesh(i, j) = 20.0    
for i = 1:n, j = 1:n
    if (i - 51)^2 + (j - 51)^2 < 50^2
        mesh[i, j]= 20.0 
    end
end

# forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 71)**2 < 10**2) mesh(i, j) =  0.0
for i = 1:n, j = 1:n
    if (i - 51)^2 + (j - 71)^2 < 10^2
        mesh[i, j]= 0.0 
    end
end

#      do iter = 0, 10000
#        if (mod(iter, 100) == 0) then        ! character animation
#          call execute_command_line('cls')  
#          print '(2g/, (52i1))', ' iteration =', iter, int(mesh(::2,::4) / 10.0)
#        end if  
#        forall (i = 1:n, j = 1:n, mask(i, j)) mesh(i, j) = (mesh(i - 1, j) + mesh(i + 1, j) + mesh(i, j - 1) + mesh(i, j + 1)) / 4
#      end do

for iter=0:1000
    if iter % 500 == 0
        # imshow(mesh[:,:]) # NG
        # imshow(mesh[0:n+1,0:n+1])  # Good
        imshow( parent(mesh)) # Good
    end
    for i = 1:n, j = 1:n
        if mask[i, j] 
            mesh[i, j] = (mesh[i - 1, j] + mesh[i + 1, j] + mesh[i, j - 1] + mesh[i, j + 1]) / 4
        end
    end
end

laplace-1.png

 計算結果は Pyplot.imshow で表示してみました。
 2次元配列 mesh が 0-based index です。Fortran 文の一部を、コメントに残しておきました。ほとんど Fortran と変わらないですね。Fortran の forall 文の 1行はよくばりすぎで、上のソースぐらいの粒度が、私は好きです。

 以下 OffsetArray の使い方を一部紹介します。

■ OffsetArrays パッケージは Pkg.add("OffsetArrays") でインストールできます。使うには using OffsetArrays です。

■ OffsetArrayを確保する方法は二つあります。

n=100
# 1
a = OffsetArray(Float64, 0:n, 0:n) 
# 2
a0 = Array(Float64, n+1, n+1)
a = OffsetArray(a0, 0:n, 0:n) 

2は、Julia標準の 1-based index配列を、0-based添字として別名を付ける例です。 型は、このように表示されます。

a0
101×101 Array{Float64,2}:

a
OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}} with indices 0:100×0:100:

■ OffsetArray は、添字を指定して読み書きする分には、通常の配列と同じように書くことができます。
■ 配列全要素を読むには、全ての範囲を添字で指定します。上の 2 の方法で確保した OffsetArray では parent(a) でもよいです。

a[0:n,0:n]
parent(a)

 同じ値を OffsetArray の全要素に代入するには fill! を使います。

fill!( a, 100.0 )
a[:,:]=100.0 # NG

■ 標準の Julia の配列操作関数の多くは、OffsetArray には対応していないので使えません。
Arrays with custom indices には、以下の注意書きがあります。
- size を使わない。代わりに indices を使う
- 1:length(A)linearindices(A) に置換する
- length(A)length(linearindices(A)) に置換する
- Array{Int}(size(B))similar(Array{Int}, indices(B)) に置換する
- @boundscheck を使わない

■ その他の表現は OffsetArray の test を見るとよいです。
https://github.com/alsam/OffsetArrays.jl/blob/master/test/runtests.jl

以上駆け足で OffsetArray の使い方を実例で見ました。0-based index 向けのパッケージは、他にもあるようです。 Julia 本体の機能として提供されるとうれしいです。 

[追記 2016/12/23 22:00]

 OffsetArray行列の一部を取り出せますね。

julia> using OffsetArrays

julia> m0=reshape( vec(1:16), (4,4))'

4×4 Array{Int64,2}:
  1   2   3   4
  5   6   7   8
  9  10  11  12
 13  14  15  16

julia> m=OffsetArray(m0,0:3,0:3)

OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}} with indices 0:3×0:3:
  1   2   3   4
  5   6   7   8
  9  10  11  12
 13  14  15  16

julia> m[1:2,1:2]

2×2 Array{Int64,2}:
  6   7
 10  11
5
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
5
4