まあ，はっきりいって，
R の敵は Python か？
じゃあ，Python の敵が Julia だったら，
Julia は R の味方じゃないですか...ニッコリ
Python なら できるんだけど，R には苦手，R じゃできないことって，どんなことがありますか？
Julia だったらできることなんじゃないですか？しかも，Python より上手に。
そんなことがあったら，知らせて下さい。
以前，Python でフィッシャーの正確確率検定を行うという記事を書いた。
今回は，もう Python は捨てようと，Julia でフィッシャーの正確確率検定を行う方法について記事を書く。
(1) http://netlib.org/toms/ から，643.gz をダウンロードする。
(2) 643.gz を展開してできる 643 というファイルを開き，1965 行目あたりからあるデータ行を削除する。
(3) 5行目（process the data file と書かれている行）の次に，以下のプログラムを挿入する。
function fisher(nrow, ncol, table)
implicit none
integer i, j, ncol, nrow
double precision emin, expect, percnt, pre, prt, table(*)
intent(in) nrow, ncol, table
double precision fisher
external fexact
expect = 0.0
percnt = 0.0
emin = 0.0
call fexact (nrow, ncol, table, nrow, expect, percnt, emin, prt,
& pre)
fisher = pre
return
end
(4) コマンドラインに以下のコマンドを入力する。
$ gfortran fisher.f -o libfisher.so -shared -fPIC
(5) エラーなくコンパイルできたら，コマンドラインに以下のコマンドを入力し，_fisher_ を含む行があるか確認する。
$ nm libfisher.so
：
000000000000665f T _fisher_
：
(6) 以下のプログラムを，ファイル名 fisher.jl でセーブする。
function fisher(tbl)
nrow, ncol = size(tbl)
ccall(("fisher_", "libfisher.so"),
Float64,
(Ref{Int32}, Ref{Int32}, Ref{Float64}),
Int32(nrow), Int32(ncol), Array{Float64, 2}(tbl))
end
(7) Julia プログラムで，以下のようにして利用する。
include("fisher.jl") # fisher.jl のパス名に注意
p = fisher([1 3 6; 4 2 3])
println("p value = $p")
(8) julia の REPL だと
$ julia
: 略
julia> include("fisher.jl")
fisher (generic function with 1 method)
julia> p = fisher([1 3 6; 4 2 3])
0.36348481240121944
julia> println("p value = $p")
p value = 0.36348481240121944
R で確認するなら
julia> using RCall
julia> R"fisher.test(matrix(c(1, 3, 6, 4, 2, 3), byrow=TRUE, ncol=3))"
RObject{VecSxp}
Fisher's Exact Test for Count Data
data: matrix(c(1, 3, 6, 4, 2, 3), byrow = TRUE, ncol = 3)
p-value = 0.3635
alternative hypothesis: two.sided
julia> R"options(digits=15)"
RObject{VecSxp}
$digits
[1] 7
julia> R"fisher.test(matrix(c(1, 3, 6, 4, 2, 3), byrow=TRUE, ncol=3))$p.value"
RObject{RealSxp}
[1] 0.363484812401219
ネット上の情報は，Julia のバージョンが低かったり，新しい書き方の FORTRAN プログラムの場合（？）だったりして，そのまま使えなかった。
試行錯誤の結果，やっとまともに動く例を示すことができた。
まずは，古典的な書式で書かれた FORTRAN プログラム（関数副プログラム，サブルーチン副プログラム）
ファイル名 test.f
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 関数副プログラム
function test1(nrow, ncol, table) ! 関数名が戻り値を表す
implicit none
integer i, j
integer nrow, ncol ! 配列のサイズ情報は引数で渡す
double precision table(*) ! 配列の受けとり方法は使い方次第（一次元配列の場合）
intent(in) nrow, ncol, table ! 入力引数
double precision test1 ! 戻り値（関数副プログラム名）の型の宣言
do 20 i=1, nrow
write (*, 99997) (table(i+(j-1)*nrow), j=1,ncol) ! 一次元配列として使用
20 continue
99997 format (1x, 10f7.1)
test1 = 0.0 ! 戻り値の定義
do 10 i = 1, nrow * ncol
test1 = test1 + table(i)
10 continue
test1 = test1 / (nrow * ncol)
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C サブルーチン副プログラム
subroutine test2(nrow, ncol, table, mean, cells) ! 戻り値は引数に指定
implicit none
integer i, j
integer nrow, ncol ! 配列のサイズ情報は引数で渡す
double precision table(nrow, *) ! 配列のサイズは最終次元のみ * でよい
double precision mean ! 戻り値の宣言
integer cells ! 戻り値の宣言
intent(in) nrow, ncol, table ! 入力引数
intent(out) mean, cells ! 戻り値の引数
mean = 0.0 ! 戻り値の設定
cells = nrow * ncol ! 戻り値の設定
do 20 i = 1, nrow ! 配列要素の並び順は Julia と同じ（列優先）
write (*, 99997) (table(i, j), j = 1, ncol) ! 二次元配列として使用
99997 format (1x, 10f7.1)
do 10 j = 1, ncol
mean = mean + table(i, j)
10 continue
20 continue
mean = mean / cells
end
次に，test.f をコンパイルし，動的ライブラリ（libtest.so）を生成する。
$ gfortran test.f -o libtest.so -shared -fPIC
できあがった動的ライブラリの中身を確認する。
$ nm libtest.so
出力される中に
0000000000000af4 T _test1_
0000000000000cd5 T _test2_
のような行がある。
_test1_，_test2_ を Julia プログラムの ccall 関数の第 1 引数として，("test1_", "libtest.so") ，("test2_", "libtest.so") として使う。
できあがった動的ライブラリを以下のような Julia プログラム test.jl から呼び出す。
# FORTRAN 関数副プログラムの例
function FuncProg(tbl)
n, m = size(tbl) # 配列サイズは引数で渡さなければならない
nrow = Int32(n) # 入力引数 Int64 を Int32 に変換
ncol = Int32(m) # 入力引数 Int64 を Int32 に変換
table = Array{Float64}(tbl) # 入力引数 Array{Int64,2}, を Array{Float64,2} に変換
# 入力引数 Int32, Int32, Array{Float64,2}
# 戻り値 Float64 FORTRAN での double precision あるいは real*8
mean = ccall(("test1_", "libtest.so"), # nm コマンドでは '_test1_' となっている
Float64, # FORTRAN 関数プログラムからの戻り値の型
(Ref{Int32}, Ref{Int32}, Ref{Float64}), # 引数の型
nrow, ncol, table) # 実引数
println("mean = $mean")
end
# FORTRAN サブルーチン副プログラムの例
function SubProg(tbl)
n, m = size(tbl) # 配列サイズは引数で渡さなければならない
nrow = Int32(n) # 入力引数 Int64 を Int32 に変換
ncol = Int32(m) # 入力引数 Int64 を Int32 に変換
table = Array{Float64}(tbl) # 入力引数 Array{Int64,2}, を Array{Float64,2} に変換
mean = [0.0] # 戻り値，スカラーのときでも必ずベクトル(値は何でもよい)
cells = [Int32(0)] # 戻り値は Int32 なので，受け取る変数も Int32 で定義
# 入力引数: Int32, Int32, Array{Float64,2} と 戻り値 Float64, Int32
# FORTRAN からの戻り値はないので Nothing を指定
ccall(("test2_", "libtest.so"),
Nothing,
(Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int32}),
nrow, ncol, table, mean, cells)
println("cells = $(cells[1]) mean = $(mean[1])") # スカラーとして取り出すときは [1] などとする
end
# FORTRAN プログラムへの引数
tbl = [1 2 3; 4 5 6]
# FORTRAN 関数副プログラム test1 を呼び出す FuncProg() の使用例
FuncProg(tbl)
# FORTRAN サブルーチン副プログラム test2 を呼び出す SubProg() の使用例
SubProg(tbl)
例えば以下のようにコマンドラインから実行する。
$ julia test.jl
結果は次のようになる。
1.0 2.0 3.0
4.0 5.0 6.0
mean = 3.5
1.0 2.0 3.0
4.0 5.0 6.0
cells = 6 mean = 3.5
