[Julia 0.5] 배열 첨자를 0부터 세기

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: htps : // 기 st. 기주 b. 이 m/그 ny 벌써 s/7437f799d132fb15에4cf0에086d5c16아9

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

계산 결과는 Pyplot.imshow로 표시해 보았습니다.
2차원 배열mesh이 0-based index입니다. Fortran 문의 일부를 주석에 남겨 둡니다. 거의 Fortran과 다르지 않네요. Fortran forall 문장의 한 줄은 너무 훌륭하고 소스가 너무 작지만 나는 좋아합니다.

이하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 는, 첨자를 지정해 읽고 쓰기하는 분에는, 통상의 배열과 같이 쓸 수가 있습니다.
■ 배열 모든 요소를 ​​읽으려면 모든 범위를 첨자로 지정합니다. 위의 두 가지 방법으로 예약한 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를 보면 좋다.
htps : // 기주 b. 이 m/아 1 m/오 f 세타라 ys. jl / b ぉ b /까지 r / st / rst sts. 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