裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

Julia の IterableTables

2021年01月31日 | ブログラミング

== IterableTables.jl

https://www.queryverse.org/IterableTables.jl/stable/

IterableTables.jl は,Julia のテーブル型データに包括的なインターフェースを定義する。

現在のところ,以下のデータソースをサポートしている。

  • DataFrames
  • DataStreams (CSV, Feather, SQLite, ODBC も含む)
  • DataTables
  • IndexedTables
  • TimeSeries
  • TypedTables
  • DifferentialEquations(あらゆる DESolution も)
  • および NamedTuple 型の要素を生成するあらゆるイテレータ

現在のところ,以下のデータシンクがサポートされている。

  • DataFrames(ModelFrame なども含む)
  • DataStreams (CSV, Feather も含む)
  • DataTables
  • IndexedTables
  • TimeSeries
  • TypedTables
  • StatsModels
  • Gadfly(currently not working)
  • VegaLite

このパッケージは Query.jl と緊密に統合されている。最後の @select 文で名前付きタプルを生成する(クエリーの結果をデータ構造に @collect しない)ような,あらゆるクエリーは自動的にイテラブルなテーブルデータソースになり,上に述べたようなデータソースは Query.jl を使って,クエリーの対象になる。

== さあ始めよう

IterableTables はJulia の異なるテーブル型の間で相互に容易に変換できる。

using IterableTables
using DataFrames

df = DataFrame(Name=["John", "Sally", "Jim"],
               Age=[34.,25.,67.],
               Children=[2,0,3])

3×3 DataFrame
 Row │ Name    Age      Children 
     │ String  Float64  Int64    
─────┼───────────────────────────
   1 │ John       34.0         2
   2 │ Sally      25.0         0
   3 │ Jim        67.0         3

using DataTables, TypedTables, IndexedTables

# DataFrame を DataTable に変換する
dt = DataTable(df)

3x3 DataTable
Name  │ Age  │ Children
──────┼──────┼─────────
John  │ 34.0 │ 2       
Sally │ 25.0 │ 0       
Jim   │ 67.0 │ 3       

# DataFrame を TypedTable に変換する
tt = Table(df)

Table with 3 columns and 3 rows:
     Name   Age   Children
   ┌──────────────────────
 1 │ John   34.0  2
 2 │ Sally  25.0  0
 3 │ Jim    67.0  3

# TypedTable を DataFrame に変換する
new_df = DataFrame(tt)

3×3 DataFrame
 Row │ Name    Age      Children 
     │ String  Float64  Int64    
─────┼───────────────────────────
   1 │ John       34.0         2
   2 │ Sally      25.0         0
   3 │ Jim        67.0         3

# TypedTable を DataTable に変換する
new_dt = DataTable(tt)

3x3 DataTable
Name  │ Age  │ Children
──────┼──────┼─────────
John  │ 34.0 │ 2       
Sally │ 25.0 │ 0       
Jim   │ 67.0 │ 3   

伝統的には DataFrame が期待されるような場合でも,どのようなデータ型も使えるようにする。

using GLM, DataFrames

# 回帰分析で TypedTable を使う
lm(@formula(Children~Age), tt)

StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Children ~ 1 + Age

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)  -0.867076    1.53833    -0.56    0.6732  -20.4134    18.6792
Age           0.0603272   0.0336492   1.79    0.3239   -0.367227   0.487881
───────────────────────────────────────────────────────────────────────────

# 回帰分析で DataTable を使う
lm(@formula(Children~Age),dt)

StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Children ~ 1 + Age

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)  -0.867076    1.53833    -0.56    0.6732  -20.4134    18.6792
Age           0.0603272   0.0336492   1.79    0.3239   -0.367227   0.487881
───────────────────────────────────────────────────────────────────────────

これらのデータソースのどれでも VegaLite でプロットできる。

using VegaLite

# TypedTable をプロットする
tt |> @vlplot(:point, x=:Age, y=:Children)

# DataTable をプロットする
dt |> @vlplot(:point, x=:Age, y=:Children)

コメント

Julia で k-means クラスター分析

2021年01月30日 | ブログラミング

まず,驚くところは,Julia の k-means 関数では,データを ncol x nrow であたえるところ。
よく,help も読まないでやると,戸惑う所の騒ぎではない。
ということで,データを transpose して与えないといけない。
これくらいのこと,ユーザに強いるなよ!!

using RDatasets
iris = dataset("datasets", "iris");
a = Matrix(iris[!, 1:2]);
using Clustering, Plots
ncluster = 3;
R = kmeans(a', 3;  maxiter=200) # データを transpose して与えること!!! a'  よりは,わかりやすく transpose(a) とする
a = assignments(R);
c = counts(R);
M = R.centers;
println(a) # 結果として,どこに分類されたか
[3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 ...
M[1,:] # 各クラスターの平均値
3-element Array{Float64,1}:
 5.800000000000001
 6.823913043478258
 5.003921568627451
図に描いてみよう
repeateach(x, n) = vec([j for i = 1:n, j in x])
colors = repeateach([:brown, :blue, :red], 50);
p1 = scatter(iris[!, 1], iris[!,2], markercolor=colors, tick_direction=:out, label="");
p1 = scatter!(M[1,:], M[2,:], markershape=:star6, markersize=10, markercolor=:cyan, label="")
display(p1)
 
 
分類結果
using FreqTables
freqtable(iris[!, 5], a)
Out[7]:
3×3 Named Array{Int64,2}
 Dim1 ╲ Dim2 │  1   2   3
────────────────────────────
setosa        │  0   0  50
versicolor    │ 38  12   0
virginica     │ 15  34   1
 R ではどんな風にやるのかな??
using RCall
R"""
a = kmeans(iris[1:2], 3)
plot(iris[,1], iris[,2], col=rep(1:3, each=50))
points(a$centers, pch=17, col=c(1,3,2), cex=2)
"""
RObject{IntSxp}
            a$cluster
iris[, 5]     1  2  3
  setosa     50  0  0
  versicolor  0 38 12
  virginica   0 15 35
コメント

Julia のデータフレームで Query

2021年01月30日 | ブログラミング
直にデータフレームを操作するより楽なのかな?
using Query
using RDatasets
mtcars = dataset("datasets", "mtcars");
names(mtcars)
12-element Array{String,1}:
 "Model"
 "MPG"
 "Cyl"
 "Disp"
 "HP"
 "DRat"
 "WT"
 "QSec"
 "VS"
 "AM"
 "Gear"
 "Carb"

@collect で,条件に合ったものをデータフレームとして返す。

q1 = @from i  in mtcars begin
    @where i.AM == 1
    @select {i.Model, i.Cyl, i.MPG, i.WT}
    @collect DataFrame
end
13 rows × 4 columns
  Model Cyl MPG WT
  String Int64 Float64 Float64
1 Mazda RX4 6 21.0 2.62
2 Mazda RX4 Wag 6 21.0 2.875
3 Datsun 710 4 22.8 2.32
4 Fiat 128 4 32.4 2.2
5 Honda Civic 4 30.4 1.615
6 Toyota Corolla 4 33.9 1.835
7 Fiat X1-9 4 27.3 1.935
8 Porsche 914-2 4 26.0 2.14
9 Lotus Europa 4 30.4 1.513
10 Ford Pantera L 8 15.8 3.17
11 Ferrari Dino 6 19.7 2.77
12 Maserati Bora 8 15.0 3.57
13 Volvo 142E 4 21.4 2.78
 

普通のデータフレームになる。

using Statistics
cor(q1[!, 3], q1[!, 4])
-0.9089147883714981
cor(Matrix(q1[:, 2:4]))
3×3 Array{Float64,2}:
  1.0       -0.825998   0.847024
 -0.825998   1.0       -0.908915
  0.847024  -0.908915   1.0

@collect しなければ,イテレータを返す。

q2 = @from i  in mtcars begin
    @where i.AM == 1
    @select {i.Model, i.WT}
end
Model WT
"Mazda RX4" 2.62
"Mazda RX4 Wag" 2.875
"Datsun 710" 2.32
"Fiat 128" 2.2
"Honda Civic" 1.615
"Toyota Corolla" 1.835
"Fiat X1-9" 1.935
"Porsche 914-2" 2.14
"Lotus Europa" 1.513
"Ford Pantera L" 3.17

... with more rows.

イテレータは for や 内包表記で使える。

for i in q2
    println("$(i.Model): weight is $(i.WT)")
end
Mazda RX4: weight is 2.62
Mazda RX4 Wag: weight is 2.875
Datsun 710: weight is 2.32
Fiat 128: weight is 2.2
Honda Civic: weight is 1.615
Toyota Corolla: weight is 1.835
Fiat X1-9: weight is 1.935
Porsche 914-2: weight is 2.14
Lotus Europa: weight is 1.513
Ford Pantera L: weight is 3.17
Ferrari Dino: weight is 2.77
Maserati Bora: weight is 3.57
Volvo 142E: weight is 2.78
[i.Model for i in q2 if i.WT > 2.5]
6-element Array{String,1}:
 "Mazda RX4"
 "Mazda RX4 Wag"
 "Ford Pantera L"
 "Ferrari Dino"
 "Maserati Bora"
 "Volvo 142E"
コメント

Julia と R のコラボレーションを jupyter notebook で記録

2021年01月30日 | ブログラミング
jupyter notebook というのは,いろいろな言語に対応している。
R の RStudio や,これまた多くの言語に対応している Atom と少し違って,スクリプトと実行結果をまとめて記録でき,成果物を LaTeX, PDF をはじめとして様々な形式で保存・閲覧できるという点で気に入っている。
たとえば,RStudio は,プログラム編集とコンソールとPlotなどのペインが分かれているので,作業しているときに煩わしい(別のやり方があるのかも知れないが)。その点,jupyter notebook はやったとおり,できたとおりのものが順に表示されていくので心地よい。あとで,不必要な部分や,追加処理を挿入するなどしてまとめて全部再度実行するのも簡単にできる。
 
それはさておき,今回は,Julia と R の環境を行き来して,いいとこ取りの分析環境を構築できるかどうか,ちょっとだけつまみ食いしてみよう。
以下では,jupiter notebook の結果を text ファイルに出力したまま(都合で html ファイルの一部分を切り取って)表示しているので不細工だが,本当は遙かに綺麗(そのまま出版できるレベル)なので,誤解なきよう。
 
まずは,Julia で,独立二標本の平均値の差の検定を行ってみる。
使うのは,ちょっと,長い名前だが UnequalVarianceTTest() だ。
 
In [1]:
using HypothesisTests
In [2]:
x = [2,1,3,2,4,3,2,3,4,5];
In [3]:
y = [3,2,3,4,5,6,3,2,3,3,3,4,3,2,6];
In [4]:
UnequalVarianceTTest(x, y)
Out[4]:
Two sample t-test (unequal variance)
------------------------------------
Population details:
    parameter of interest:   Mean difference
    value under h_0:         0
    point estimate:          -0.566667
    95% confidence interval: (-1.6209, 0.4876)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.2760

Details:
    number of observations:   [10,15]
    t-statistic:              -1.1192102478745312
    degrees of freedom:       20.56776717709893
    empirical standard error: 0.506309397848002
次に,Julia から RCall を介して,R の t.test() で同じデータを分析してみる
 
In [5]:
using RCall
長いスクリプトは,R""" と """ で挟んで,何行でも書ける。
$x, $y というのは Julia の x, y ということを表す。
結果を明示的に print() で書かないと,最後の出力しか表示されないので注意。

In [6]:
R"""
options(digits=15)
a = t.test($x, $y)
print(a)
print(a$statistic)
print(a$p.value)
"""
 
	Welch Two Sample t-test

data:  `#JL`$x and `#JL`$y
t = -1.119210247875, df = 20.5677671771, p-value = 0.275951369175
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.620943659272958  0.487610325939624
sample estimates:
       mean of x        mean of y 
2.90000000000000 3.46666666666667 

                t 
-1.11921024787453 
[1] 0.275951369174595
Out[6]:
RObject{RealSxp}
[1] 0.275951369174595
In [7]:
R"a"
Out[7]:
RObject{VecSxp}

	Welch Two Sample t-test

data:  `#JL`$x and `#JL`$y
t = -1.119210247875, df = 20.5677671771, p-value = 0.275951369175
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.620943659272958  0.487610325939624
sample estimates:
       mean of x        mean of y 
2.90000000000000 3.46666666666667 

図を描くと,図もその場に出力されるので快適。
ただ,希望しない出力もしてくれるのは大きなお世話ともいえるが,普通は気にしない分析結果なので,統計解析もしっかりやらなければならないよという戒めにもなるのかな?

In [8]:
R"""
group = rep(factor(c("x", "y")), c(length($x),length($y)))
value = c($x, $y)
a = boxplot(value ~ group)
"""
Out[8]:
RObject{VecSxp}
$stats
     [,1] [,2]
[1,]    1    2
[2,]    2    3
[3,]    3    3
[4,]    4    4
[5,]    5    5
attr(,"class")
        x 
"integer" 

$n
[1] 10 15

$conf
                 [,1]             [,2]
[1,] 2.00072025938679 2.59204575419949
[2,] 3.99927974061321 3.40795424580051

$out
[1] 6 6

$group
[1] 2 2

$names
[1] "x" "y"

コメント

汎用性のあるプログラム(関数)を書こうよ

2021年01月29日 | ブログラミング

毎度,おいしいテーマを提供いただき,ありがとうございます。

「全て異なる」という条件をPythonで実験する

ということですが,

要素数を無制限にした,以下のような回答があるかと思いきや。

a = [1,2,3,4,5]
all([a[i] != a[j] for i in range(len(a)-1) for j in range(i+1,len(a))]) # True

a = [1,2,1,4,5]
all([a[i] != a[j] for i in range(len(a)-1) for j in range(i+1,len(a))]) # False

残念でした。

Julia だと,

allequal(a) = all([a[i] != a[j] for j = 1:length(a), i = 1:length(a) if i < j])

という関数を作ればいつでも使えます。

a = [1,2,3,4,5]
allequal(a) # false
a = [1,2,1,4,5]
allequal(a) # true

いやいやそんな冗長な!!(^_^;)

length(Set(a)) == length(a)

 で,いいじゃん!!チャンチャン!

コメント

Julia で R を使う(その2)

2021年01月29日 | ブログラミング

== サポートされている変換

RCall は最も基本的な Julia のデータ型や人気のある統計パッケージ例えば DataFrames, CategoricalArrays や AxisArrays との相互変換を支援する。

== Julia の基本データ型

julia> # Julia -> R
julia> a = robject(1)
RObject{IntSxp}
[1] 1

julia> # R -> Julia
julia> rcopy(a)
1

julia> # Julia -> R
julia> a = robject([1.0, 2.0])
RObject{RealSxp}
[1] 1 2

julia> # R -> Julia
julia> rcopy(a)
2-element Array{Float64,1}:
 1.0
 2.0

== 辞書型 Dictionaries

julia> # Julia -> R
julia> d = Dict(:a => 1, :b => [4, 5, 3])
julia> r = robject(d)
RObject{VecSxp}
$a
[1] 1

$b
[1] 4 5 3

julia> # R -> Julia
julia> rcopy(r)
OrderedCollections.OrderedDict{Symbol,Any} with 2 entries:
  :a => 1
  :b => [4, 5, 3]

== 日付 Date

julia> using Dates
julia> # Julia -> R
julia> d = Date(2012, 12, 12)
julia> r = robject(d)
RObject{RealSxp}
[1] "2012-12-12"

julia> # R -> Julia
julia> rcopy(r)
2012-12-12

== 日付と時間 DateTime

julia> # julia -> R
julia> d = DateTime(2012, 12, 12, 12, 12, 12)
julia> r = robject(d)
RObject{RealSxp}
[1] "2012-12-12 12:12:12 UTC"

julia> # R -> Julia
julia> rcopy(r)
2012-12-12T12:12:12

== データフレーム DataFrames

julia> using DataFrames
julia> d = DataFrame([[1.0, 4.5, 7.0]], [:x])
julia> # Julia -> R
julia> r = robject(d)
RObject{VecSxp}
    x
1 1.0
2 4.5
3 7.0

julia> # R -> Julia
julia> rcopy(r)
3×1 DataFrame
 Row  │ x       
      │ Float64 
──────┼─────────
   1  │     1.0
   2  │     4.5
   3  │     7.0

デフォルトでは,R のデータフレームの列名(変数名)は foo.bar が foo_bar になるような正規化が行われる。

julia> rcopy(R"data.frame(a.b = 1:3)")
3×1 DataFrame
 Row  │ a_b   
      │ Int64 
──────┼───────
   1  │     1
   2  │     2
   3  │     3

正規化を避けるためには,normalizenames オプションを使う。

julia> rcopy(R"data.frame(a.b = 1:10)"; normalizenames = false)
10×1 DataFrame
 Row  │ a.b   
      │ Int64 
──────┼───────
   1  │     1
   2  │     2
   3  │     3
   4  │     4
   5  │     5
   6  │     6
   7  │     7
   8  │     8
   9  │     9
  10  │    10
  
== 名前付き配列 AxisArrays

julia> using AxisArrays
julia> # Julia -> R
julia> aa = AxisArray([1,2,3], Axis{:id}(["a", "b", "c"]))
julia> r = robject(aa)
RObject{IntSxp}
id
a b c 
1 2 3 

julia> # R -> Julia
julia> rcopy(AxisArray, r)
1-dimensional AxisArray{Int64,1,...} with axes:
    :id, ["a", "b", "c"]
And data, a 3-element Array{Int64,1}:
 1
 2
 3

julia> bb = AxisArray([1 2 3;4 5 6;7 8 9], Axis{:id}(["a", "b", "c"]), Axis{:col}(["A1", "A2", "B"]))
julia> s = robject(bb)
RObject{IntSxp}
   col
id  A1 A2 B
  a  1  2 3
  b  4  5 6
  c  7  8 9

julia> # R -> Julia
julia> rcopy(AxisArray, s)
2-dimensional AxisArray{Int64,2,...} with axes:
    :id, ["a", "b", "c"]
    :col, ["A1", "A2", "B"]
And data, a 3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

== 特別な変換

RCall はR と Julia のオブジェクトの暗黙的な変換を rcopy と robject で行うのをサポートするための API を提供する。
このアイデアを説明するために,以下のような Julia のデータ型を考える。

julia> mutable struct Foo
           x::Float64
           y::String
       end

julia> foo = Foo(1.0, "hello")
julia> bar = R"""
           bar <- list(x = 1, y = "hello")
           class(bar) <- "Bar"
           bar
       """
RObject{VecSxp}
$x
[1] 1

$y
[1] "hello"

attr(,"class")
[1] "Bar"

== R から Julia へ

rcopy 関数と rcopytype 関数はこの方向の変換を受け持つ。

julia> import RCall.rcopy

julia> function rcopy(::Type{Foo}, s::Ptr{VecSxp})
           Foo(rcopy(Float64, s[:x]), rcopy(String, s[:y]))
       end
rcopy (generic function with 92 methods)

convert 関数があると対応する rcopy 関数を起動する。

julia> rcopy(Foo, bar)
julia> convert(Foo, bar) # calls `rcopy`
Foo(1.0, "hello")

rcopy(var)が自動的に変換するのを許可するためには,R のクラス Barが登録されていなければならない。

julia> import RCall: RClass, rcopytype
       
julia> rcopytype(::Type{RClass{:Bar}}, s::Ptr{VecSxp}) = Foo
julia> foo2 = rcopy(bar)
Foo(1.0, "hello")

== Julia から R へ

Julia から R への変換を許可するためには RCall.sexp 関数は上書き(再定義)されなければならない。
sexp 関数は Julia オブジェクトを受け付け,SEXP オブジェクト( [Sxp] へのポインタ)を返す。
まず最初に,Julia 型の Foo から R のクラス Bar への明示的な変換法を定義する。

julia> import RCall: sexp, protect, unprotect, setclass!, RClass
 
julia> function sexp(::Type{RClass{:Bar}}, f::Foo)
           r = protect(sexp(Dict(:x => f.x, :y => f.y)))
           setclass!(r, sexp("Bar"))
           unprotect(1)
           r
       end
       
julia> bar = robject(:Bar, foo)
RObject{VecSxp}
$y
[1] "hello"

$x
[1] 1

attr(,"class")
[1] "Bar"

注:SEXP がガーベージコレクションされないように保護するために RCall.protect と RCall.unprotect を使わなければならない。
robject(foo) を経由するデフォルトの変換を登録するために,sexpclass を定義する必要がある。

julia> import RCall.sexpclass
 
julia> sexpclass(f::Foo) = RClass{:Bar}
julia> bar = robject(foo)
RObject{VecSxp}
$y
[1] "hello"

$x
[1] 1

attr(,"class")
[1] "Bar"

== @rput と @rget をシームレスに使う

julia> foo2.x = 2.0
julia> @rput foo2
julia> R"""
       foo2["x"]
       """
RObject{VecSxp}
$x
[1] 2

julia> R"""
       foo2["x"] = 3.0
julia> """
julia> @rget foo2
julia> foo2.x
3.0

== 入れ子にされた変換

julia> l = R"list(foo2 = foo2, bar = $bar)"
RObject{VecSxp}
$foo2
$y
[1] "hello"

$x
[1] 3

attr(,"class")
[1] "Bar"

$bar
$y
[1] "hello"

$x
[1] 1

attr(,"class")
[1] "Bar"

julia> rcopy(l)
OrderedCollections.OrderedDict{Symbol,Any} with 2 entries:
  :foo2 => Foo(3.0, "hello")
  :bar  => Foo(1.0, "hello")

== イベントループ

IJulia ではない対話セッションでは,R はデフォルトで新たな plot 表示ウインドウを開く。
プロットのサイズ変更などを対話的に行うために,R イベントループが自動的に開始される。
以下のようにすれば,手動で R のイベントループを開始できる。

julia> RCall.rgui_start()
true

これにより,プロットが変更されたかどうかをチェックするために R を頻繁に呼出すことになる。
終了するためには,以下のようにする。

julia> RCall.rgui_stop()
true

 

コメント

Julia で R を使う(その1)

2021年01月29日 | ブログラミング

== RCall.jl

参照 URL https://juliainterop.github.io/RCall.jl/stable/

RCall は,Julia から R 環境(R パッケージ)を利用したり,R 環境での作業結果を Julia へ取り込んだりすることができる。

== インストール

まだ RCall をインストールしていない場合は,以下の 2 行を実行。

julia> using Pkg
julia> Pkg.add("RCall")

== さあ始めよう

まず,最初に RCall パッケージをロードしよう。
他のパッケージと同じで,あるパッケージを使いたいときにはそのセッションの最初に「using パッケージ名」と指定する。

julia> using RCall

== RCall の使い方

Julia と R の相互作用のために,RCall は複数の方法を提供する。

  • R の REPL モード
  • @rput と @get マクロ
  • R"文字列" マクロ
  • RCall API: reval, rcall, rcopy, robject など
  • @rlibrary と @rimport マクロ

== R の REPL モード

Julia のプロンプトが出ているときに $ を押すと,R の REPL モードに移行する(プロンプトが julia> から R> に変わる)。
R から Julia に戻るときは delete キー(Windows の場合は backspaceキー)を」押す(プロンプトがR> から julia> に戻る)。
R の REPL モードで Julia の変数(オブジエクト)を引用するときには「$変数名」,「$]式」のようにする。
julia の println() で,文字列中に $変数名 と書けばその値で置き換えられるのと同じような使い方である。

julia> foo = 1
1

$ キーを押す

R> x <- $foo

R> x
[1] 1

R> y = $(rand(10))

R> sum(y)
[1] 5.234343

delete キーを押す

julia>

== @rput マクロと @rget マクロ

これらのマクロは,Julia と R の間で変数を転送する。
コピーされた変数は元の変数と同じ名前になる。

julia> z = 1
1

julia> @rput z
1

R> z
[1] 1

R> r = 2

julia> @rget r
2.0

julia> r 2.0

一行で複数の変数をやり取りすることもできる。

julia> foo = 2
2

julia> bar = 4
4

julia> @rput foo bar
4

R> foo + bar
[1] 6

== @R_str 文字列マクロ

RCall の別の使い方は R"文字列" マクロである。
これは特にスクリプトファイルで有用である。

julia> R"rnorm(10)" # プロンプトは julia> の状態で使う
RObject{RealSxp}
 [1]  2.12249461  0.09807257 -2.28437327 ...
 [7] -0.49853088  1.16451855 -0.41658982 ...

このような使い方をすると,文字列を R の式とみなし,R での評価結果を Julia のラッパーでくるんだ RObject として返す。
R"文字列" マクロは変数代入(\latexmath:[$変数名)をサポートする。
aa\\$]bb のような R の式としては不適切な場合にさえも使うことができる。

julia> x = randn(10)
10-element Array{Float64,1}:
 -2.8182219563281503
  0.8334083608216333
 -1.2866572484420662
 -0.719998093988448
  1.3888582310416397
  0.45538366532280794
  0.26108314692800977
 -0.2748229904314617
 -0.9151614036341179
  1.1915882754899525

julia> R"t.test($x)"
RObject{VecSxp}

 One Sample t-test

data:  `#JL`$x
t = -0.46172, df = 9, p-value = 0.6552
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -1.1117626  0.7348546
sample estimates:
mean of x 
-0.188454 

Julia で評価される前の式を R に渡すこともできる。
このような場合には ( ) でくくる必要がある。
なお,以下の例では文字列をくくるのに二重引用符 " をつかっているが,その中の文字列中に出てくる文字列は一重引用符 ’を使う例を示している。
一重引用符の代わりに「エスケープ付の二重引用符 \"」を使ってもよい。
外側の引用符は二重引用符でなければならない(一重引用符を使うとエラーになる)。

julia> R"optim(0, $(x -> x-cos(x)), method='BFGS')"
RObject{VecSxp}
$par
[1] -1.56343

$value
[1] -1.570796

$counts
function gradient 
      14       13 

$convergence
[1] 0

$message
NULL

複数行にわたるような文字列スクリプトは,(Python と同じ)三重引用符(二重引用符を3個連続)を使う。

julia> y = 1
1

julia> R"""
       f <- function(x, y) x + y
       ret <- f(1, $y)
       """

RObject{RealSxp}
[1] 2

== RCall API

reval 関数は,入力された文字列を R コードとして評価する。
評価結果は RObject として返される。

julia> jmtcars = reval("mtcars");
julia> typeof(jmtcars)
RObject{VecSxp}

julia> names(jmtcars)
11-element Array{Symbol,1}:
 :mpg
 :cyl
 :disp
 :hp
 :drat
 :wt
 :qsec
 :vs
 :am
 :gear
 :carb

julia> jmtcars[:mpg]
RObject{RealSxp}
 [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 ...
[16] 10.4 14.7 32.4 30.4 33.9 21.5 15.5 15.2 ...
[31] 15.0 21.4

julia> df = rcopy(jmtcars); # データフレームに変換
julia> first(df, 5)
5×11 DataFrame
 Row │ mpg      cyl      disp     hp       drat ...
     │ Float64  Float64  Float64  Float64  Float64 ...
─────┼──────────────────────────────────────────────────
   1 │    21.0      6.0    160.0    110.0     3.9 ...
   2 │    21.0      6.0    160.0    110.0     3.9 ...
   3 │    22.8      4.0    108.0     93.0     3.85 ...
   4 │    21.4      6.0    258.0    110.0     3.08 ...
   5 │    18.7      8.0    360.0    175.0     3.15 ...

rcall 関数は関数呼び出しのために使用される。

jmtcars の次元サイズを求め RObject に変換する
julia> rcall(:dim, jmtcars)
RObject{IntSxp}
[1] 32 11

引数は暗黙的に評価され RObject に変換される。

ベクトルの和を求め RObject に変換する
julia> rcall(:sum, Float64[1.0, 4.0, 6.0])
RObject{RealSxp}
[1] 11

rcopy 関数はRObject を Julia のオブジェクトに変換する。
豊富な発見的手法に基づいて,もっとも適切な Julia の型に変換する。

julia> rcopy(R"c(1)")
1.0

julia> rcopy(R"c(1, 2)")
2-element Array{Float64,1}:
 1.0
 2.0

julia> rcopy(R"list(1, 'zz')")
2-element Array{Any,1}:
 1.0
 "zz"

julia> rcopy(R"list(a = 1, b= 'zz')")
OrderedCollections.OrderedDict{Symbol,Any} with 2 entries:
  :a => 1.0
  :b => "zz"

特定の型変換をするために,第 1 引数として型を指定することができる。

julia> rcopy(Array{Int}, R"c(1, 2)")
2-element Array{Int64,1}:
 1
 2

robject 関数はあらゆる Julia オブジェクトを RObject に変換する。

julia> robject(1)
RObject{IntSxp}
[1] 1

julia> robject(Dict(:a => 1, :b => 2))
RObject{VecSxp}
$a
[1] 1

$b
[1] 2

== @rlibrary マクロと @import マクロ

このマクロは,全てのエクスポートされた R パッケージ中の関数およびオブジェクトを現在のモジュールにロードする。

julia> @rlibrary boot
julia> city = rcopy(R"boot::city")  # get some data
10×2 DataFrame
│ Row  │ u      │ x      │
│      │ Float64│ Float64│
├──────┼────────┼────────┤
│ 1    │ 138.0  │ 143.0  │
│ 2    │ 93.0   │ 104.0  │
│ 3    │ 61.0   │ 69.0   │
│ 4    │ 179.0  │ 260.0  │
│ 5    │ 48.0   │ 75.0   │
│ 6    │ 37.0   │ 63.0   │
│ 7    │ 29.0   │ 50.0   │
│ 8    │ 23.0   │ 48.0   │
│ 9    │ 30.0   │ 111.0  │
│ 10   │ 2.0    │ 50.0   │

julia> ratio(d, w) = sum(d[!, :x] .* w)/sum(d[!, :u] .* w)
ratio (generic function with 1 method)

julia> b = boot(city, ratio, R = 100, stype = "w");
julia> rcall(:summary, b[:t])
RObject{StrSxp}
       V1       
 Min.   :1.274  
 1st Qu.:1.449  
 Median :1.543  
 Mean   :1.572  
 3rd Qu.:1.638  
 Max.   :2.178  

当然ではあるが,データは R と Julia の間で何回もコピーされるので,余りにも非効率である。
効率性の面からは,R"文字列" マクロを使うことをお勧めする。
いくつかの R 関数はドットを含むキーワード引数を持っている。
RCall はこれらのキーワードをエスケープするために,文字列マクロ「var"ドットを含む引き数名"」を提供する。

julia> @rimport base as rbase
julia> rbase.sum([1, 2, 3], var"rm.na" = true)
RObject{IntSxp}
[1] 7

コメント

MATLAB プログラムを Julia に書き換え

2021年01月29日 | ブログラミング

細かいところは最終的には頭を使った手作業になるが,エディタレベルでの作業についてまとめておく。

  • 文末が ; であるが,Julia では問題ないのでそのままにしておく。
  • if 文,while 文 の論理式は ( )  で囲まれているが,Julia では問題ないのでそのままにしておく。
  • 関数が function [戻り値リスト] = 関数名(入力引数)になっているが,"[戻り値リスト] =" を取り除く。
  • 注釈行は % で始まる % を # に全置換する。
  • プログラミングスタイルにもよるが,関数名と ( の間の空白は許されないので "関数名(" のようにする。
  • persistent 変数は必要ないので,persistent 文をコメントアウトする。
  • 文字列は ' ではなく" でくくるので,全置換する(' のままのものもあるので,確認しつつ置換する)。
  • 論理式の否定演算子が ~ なので,! に全置換する。
  • 出力文が fprintf() なので,@printf() に置換する。  "fprintf ( 1, "を "@printf(" に全置換する。先頭に "using Printf"を加えておく。println() に置換してもよいが format 付の場合には @printf() が妥当。
  • zeros(), ones() は,ベクトルであっても,二次元目が 1 になっているので,zeros(n, 1) のようなものを zeros(1),さらに型を考えて zeros(Int64, n), ones(Float64, n) のように直す。
  • ベクトル,配列は V(i),A(i, j) のようになっているので,V[i], A[i, j] にする。全置換と手作業の組み合わせ。
  • 継続行は ... で終わるが,これは単に ... を消去するだけでよい。
  • 論理値を 0/1 で表すので,false/true に変換。
  • 一様乱数の使用は [ r, seed ] = r8_uniform_01 ( seed ) のようになっているので,r = rand() のようにする。
  • 戻り値がちゃんと返されるように,必要なら return 文で指定する。
  • この辺りまで来たら,REPL で関数全体をコンパイルし,エラーメッセージに従ってプログラムを修正する。
  • 別の記事にも書いたが,Julia での for ループ,while ループに変数スコープがあるために,"変数名 not defined" に対応する。
コメント

Julia の変数スコープ

2021年01月28日 | ブログラミング

わかりにくいかもしれないので,簡単な例を挙げておこう。

test99.jl という名前のファイルに以下の内容が入っているとする

for i = 1:10
    s = i == 1
    s = s + i
end
println(s)

コマンドラインで

$ julia test99.jl

とすると,

ERROR: LoadError: UndefVarError: s not defined
Stacktrace:
 [1] top-level scope at /Users/foo/Desktop/Julia/test99.jl:5
 [2] include(::Function, ::Module, ::String) at ./Base.jl:380
 [3] include(::Module, ::String) at ./Base.jl:368
 [4] exec_options(::Base.JLOptions) at ./client.jl:296
 [5] _start() at ./client.jl:506
in expression starting at /Users/foo/Desktop/Julia/test99.jl:5

となる。

REPL で,一番最初に

include("test99.jl")

を実行すると,同じエラーが出る(少し違うが)

LoadError: UndefVarError: s not defined
in expression starting at /Users/foo/Desktop/Julia/test99.jl:5

Stacktrace:
 [1] top-level scope at /Users/foo/Desktop/Julia/test99.jl:5
 [2] include(::String) at ./client.jl:457
 [3] top-level scope at In[2]:1
 [4] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

この後例えば

s = 0

という行を実行させた後,再度

include("test99.jl")

を実行すると

0
┌ Warning: Assignment to `s` in soft scope is ambiguous because a global variable by the same name exists: `s` will be treated as a new local. Disambiguate by using `local s` to suppress this warning or `global s` to assign to the existing global variable.
└ @ nothing /Users/foo/Desktop/Julia/test99.jl:2

となる。

さて,普通のプログラミング言語だとどうなっていたっけ?と思い,まずは R でやってみる。

Julia でのプログラムを以下のように書き換えて test99.R に保存する。

for (i in 1:10) {
    s = i == 1
    s = s + i
}
print(s)

コマンドラインからこれを Rscript で実行すると

$ Rscript test99.R
[1] 10

なんのエラーもなく,期待される答え 10  を表示する。

R か RCommander で

source("test99.R")

としても,同じく 10 を表示する。エラーはない。

次に Python でやってみる。以下のプログラムを test99.py として保存する。

for i in range(1, 11):
    s = i == 1
    s = s + i
print(s)

コマンドラインからは

$ python3 test99.py
10

となるし,REPL でもエラーなく,同じ答え 10 が表示される。

 

我々としては,R や Python の振る舞いの方が普通に思われるが,確かに for ループの中で最初に s にアクセスするときは,「s なんかないよ」というのが理にかなっているのかも知れないなあ。

さて,悩んでいてもどうにもならないので,対策を考える。

for の内と外で s といっても違いがあるということで,

s = 999
for i = 1:10
    s = i == 1
    s = s + i
end
println(s)

としただけでは,まだ期待するような結果にはならない。

999
┌ Warning: Assignment to `s` in soft scope is ambiguous because a global variable by the same name exists: `s` will be treated as a new local. Disambiguate by using `local s` to suppress this warning or `global s` to assign to the existing global variable.
└ @ nothing /Users/aoki/Desktop/Julia/test99.jl:3

ところで,関数もまた変数スコープを持つので,このプログラムを関数定義でくるんでやる。

function func()
    s = 999
    for i = 1:10
        s = i == 1
        s = s + i
    end
    println(s)
end

これではじめて,for の内と外の s が同じものになる。その後で,関数を読んでやれば,期待する答え 10 が表示される。

s = 999初期化は不要だが定義が必要としたものだ。

func()
10

 

 

コメント

Julia の変数スコープに注意

2021年01月28日 | ブログラミング

R の fisher.test() において,hybrid=true を選択した場合に使用される rcont2() を Julia に移植した。

Julia は,for ブロック,while ブロックが 新たな変数スコープを持つので,注意が必要な局面がある。

普通は,エラーメッセージで気づくのだが,REPL 環境でやっているとおかしなことになってしまい,ドツボに嵌まることもある。

この関数でも,3 箇所(3 変数)で注意が必要。初期化は必要ない(どのような数値を代入しても無関係)だが,そこでその変数を定義しないと「そんな変数ないよ」エラーになる。試行錯誤・テストのために PEPL で適当な値を設定してしまっていれば,このエラーが出なくなってします。そこで,ドツボにはまる(for や while の外に出たとたんに変数の値が代わってしまうという現象が起きる)。

もうひとつは,Julia らしい,面白いプログラミング技法

論理式 && 論理式が真の場合に実行される文

論理式 || 論理式が偽の場合に実行される文

if a < 5
     break
end

なんてのが,一行で書ける(わかりにくくなってしまうかもしれないが,慣れると便利)

using SpecialFunctions

function rcont2(nrowt, ncolt)
  #=============================================================================
  RCONT2 constructs a random two-way contingency table with given sums.
  Licensing:
    This code is distributed under the GNU LGPL license.
  Modified:
    28 Janualy 2021
  Author:
    Original FORTRAN77 version by WM Patefield.
  Reference:
    WM Patefield,
    Algorithm AS 159:
    An Efficient Method of Generating RXC Tables with
    Given Row and Column Totals,
    Applied Statistics,
    Volume 30, Number 1, 1981, pages 91-97.
  Parameters:
    Input, Array{Int64,1}nrowt(nrow), Array{Int64}ncolt(ncol),
      the row and column sum vectors.
    Output, Array{Int64,2}matrix(nrow, ncol), the matrix.
  =============================================================================#

  MAXTOT = 1000000
  nrow = length(nrowt)
  ncol = length(ncolt)
  ntotal = sum(ncolt)
  nrow > 1 || error("The length of the row sum vector must be greater than 1.")
  ncol > 1 || error("The length of the column sum vector must be greater than 1.")
  all(nrowt .> 0) || error("An entry in the row sum vector must be greater than 0.")
  any(ncolt .> 0) || error("An entry in the column sum vector must be greater than 0..")
  ntotal == sum(nrowt) || error("The row and column sum vectors must have the same sum.")
  MAXTOT >= ntotal || error("The sum of the column sum vector entries is too large.")

  # Calculate log-factorials.
  fact = logfactorial.(0:ntotal)
  # Construct a random matrix.
  matrix = zeros(Int64, nrow, ncol)
  jwork = zeros(Int64, ncol)
  jwork[1:ncol-1] = ncolt[1:ncol-1]
  jc = ntotal
  ib = 99999 # 初期化は不要だが,定義は必要
  for l = 1 : nrow - 1
    nrowtl = nrowt[l]
    ia = nrowtl
    ic = jc
    jc -= nrowtl
    for m = 1 : ncol - 1
      id = jwork[m]
      ie = ic
      ic -= id
      ib = ie - ia
      ii = ib - id
      # Test for zero entries in matrix.
      if ie == 0
        ia = 0
        matrix[l, m:ncol] = 0
        break
      end
      # Generate a pseudo-random number.
      r = rand(1)[1]
      # Compute the conditional expected value of MATRIX(L,M).
      done1 = false
      nlm = 99999 # 初期化は不要だが,定義は必要
      while true
        nlm = floor(Int64, ia * id / ie + 0.5)
        iap = ia + 1
        idp = id + 1
        igp = idp - nlm
        ihp = iap - nlm
        nlmp = nlm + 1
        iip = ii + nlmp
        x = exp(fact[iap] + fact[ib+1] + fact[ic+1] + fact[idp] -
          fact[ie+1] - fact[nlmp] - fact[igp] - fact[ihp] - fact[iip])
        r <= x && break
        sumprb = x
        y = x
        nll = nlm
        lsp = false
        lsm = false
        # Increment entry in row L, column M.
        done2 =false # 初期化は不要だが,定義は必要
        while ! lsp
          j = (id - nlm) * (ia - nlm)
          if j == 0
            lsp = true
          else
            nlm += 1
            x *= j / (nlm * (ii + nlm))
            sumprb = sumprb + x
            if r <= sumprb
              done1 = true
              break
            end
          end
          done2 = false
          while !lsm
            # Decrement the entry in row L, column M.
            j = nll * (ii + nll)
            if j == 0
              lsm = true
              break
            end
            nll -= 1
            y = y * j / ((id - nll) * (ia - nll))
            sumprb = sumprb + y
            if r <= sumprb
              nlm = nll
              done2 = true
              break
            end
            !lsp && break
          end
          done2 && break
        end
        (done1 || done2) && break
        r = sumprb * rand(1)[1]
      end
      matrix[l, m] = nlm
      ia -= nlm
      jwork[m] -= nlm
    end
    matrix[l, ncol] = ia
  end
  # Compute the last row.
  matrix[nrow, 1:ncol-1] = jwork[1:ncol-1]
  matrix[nrow, ncol] = ib - matrix[nrow, ncol-1]
  return matrix
end

# 使い方

rcont2([1,2,3], [2,4])

コメント

みんな大好き,Julia の速度自慢

2021年01月25日 | ブログラミング

つい一月半位前に

みんな大好き,R と Python の速度比較

というのを書いた。そのときには,Julia なんて眼中になかったのだ。

そこで,Julia で同じことをやるとどんだけ〜〜速いかやってみた。

 

結論

Julia 速い。ただし,初回は若干時間が掛かる。

Python は ピアソンの積率相関係数行列,分散・共分散行列,重回帰分析で健闘している。

 

using CSV, DataFrames

test.csv は 100000行, 1500列のデータ

データ読み込み

@time df = CSV.read("test.csv", DataFrame)
#  23.113033 seconds (24.45 M allocations: 2.275 GiB, 3.44% gc time)
# R では,変数型を指定して 33.607, Pythonでは 28.793 だったので,かなり速いといえよう
# 以下では 単位は秒

単変量分析

@time a = mapcols(sum, df)  # 0.0777, R: 0.739, Pyhton: 5.416
using Statistics
@time a = mapcols(mean, df) # 0.0760, R: 0.662, Python:  5.352
@time a = mapcols(var, df)  # 0.123,  R: 0.552, Python:  8.466
@time a = mapcols(std, df)  # 0.125,  R: 0.561, Python: 11.574

2変量分析

ピアソンの積率相関係数
@time a = cor(df[:, 1], df[:, 2]) # 0.00106, R: 0.011, Python: 0.004

多変量回帰

ピアソンの積率相関係数行列
@time a = cor(Array{Float64,2}(df)) # 4.331, R: 100.509, Python: 360.236

スピアマンの順位相関係数行列
using StatsBase
function spearmancor(x::Array{Float64,2})
    rank = similar(x, Float64)
    for i = 1:size(x, 2)
        rank[:, i] = tiedrank(x[:, i])
    end
    return cor(rank)
end
@time a = spearmancor(Array{Float64,2}(df)) # 23.354, R: 140.892, Python 4238.496

分散・共分散行列
@time a = cov(Array{Float64,2}(df));  # 4.469, R: 106.829, Python: 4.987

重回帰分析
using GLM
@time a = lm(@formula(X1 ~ X2+X3+X4+X5+X6+X7+X8+X9+X10), df) # 0.0255, R: 0.163, Python 0.031

コメント

簡単なプログラミング問題を Julia で解く

2021年01月25日 | ブログラミング

以下では Python の解答例が示されているので,Julia でやったらどうなるか書いてみた。

10 Algorithms To Solve Before your Python Coding Interview

関数として定義したが,一般的に Python よりは短く簡潔に書ける。

1. Reverse Integer

Given an integer, return the integer with reversed digits. 
Note: The integer could be either positive or negative.

func1(x) = (x < 0 ? "-" : "") * reverse(string(abs(x)))

println(func1(-123))
println(func1(123))

-321
321

2. Average Words Length

For a given sentence, return the average word length. Note: Remember to
remove punctuation first.

たかが mean() を使うだけで using Statistics するのは,大げさだなあ。

using Statistics
func2(sentence) = mean(length.(split(replace(sentence, r"[!?':,,.]" => ""))))

sentence1 = "Hi all, my name is Tom... I am originally from Australia."
sentence2 = "I need to work very hard to learn more about algorithms in Python!"
println(func2(sentence1))
println(func2(sentence2))

3.8181818181818183
4.076923076923077

3. Add Strings

Given two non-negative integers num1 and num2 represented as string,
return the sum of num1 and num2. You must not use any built-in
BigInteger library or convert the inputs to integer directly. Notes:
Both num1 and num2 contains only digits 0-9. Both num1 and num2 does not
contain any leading zero

func3(a, b) = parse(Int, a) + parse(Int, b)

num1 = "364"
num2 = "1836"
println(func3(num1, num2))

2200

4. First Unique Character

Given a string, find the first non-repeating character in it and return
its index. If it doesn’t exist, return -1. Note: all the input strings
are already lowercase.

function func4(s)
    t = split(s, "")
    for i in 1:length(t)
        if count(t .t[i]) 1
            return i
        end
    end
    return -1
end

println(func4("alphabet")) # Julia では,文字位置は1から数える
println(func4("barbados"))
println(func4("ppaap"))

2
3
-1

5. Valid Palindrome

Given a non-empty string s, you may delete at most one character. Judge
whether you can make it a palindrome. The string will only contain
lowercase characters a-z.

function func5(s)
    len = length(s)
    for i = 1:len
        t = s[1:i-1] * s[i+1:len]
        if t reverse(t)
            return true
        end
    end
    return false
end

s = "radkar"
println(func5(s))
println(func5("abcde"))

true
false

6. Monotonic Array

Given an array of integers, determine whether the array is monotonic or
not.

func6(x) = x sort(x) || x sort(x, rev=true)

println(func6([6, 5, 4, 4]))
println(func6([1, 1, 1, 3, 3, 4, 3, 2, 4, 2]))
println(func6([1, 1, 2, 3, 7]))

true
false
true

7. Move Zeroes

Given an array nums, write a function to move all zeroes to the end of
it while maintaining the relative order of the non-zero elements.

発想の転換をすれば,解は簡単に求まる。

function func7(x)
    a = [i for i in x if i != 0]
    append!(a, repeat([0], count(x .0)))
    return a
end

array1 = [0,1,0,3,12]
println(func7(array1))
array2 = [1,7,0,0,8,0,10,12,0,4]
println(func7(array2))

[1, 3, 12, 0, 0]
[1, 7, 8, 10, 12, 4, 0, 0, 0, 0]

8. Fill The Blanks

Given an array containing None values fill in the None values with most
recent non None value in the array.

Julia に None はないので,NaN に置き換えて考える。

function func8(x)
    for i in 2:length(x)
        if isnan(x[i])
            x[i] = x[i-1]
        end
    end
    x
end

array1 = [1,NaN,2,3,NaN,NaN,5,NaN];
println(func8(array1))

[1.0, 1.0, 2.0, 3.0, 3.0, 3.0, 5.0, 5.0]

9. Matched & Mismatched Words

Given two sentences, return an array that has the words that appear in
one sentence and not the other and an array with the words in common.

function func9(a, b)
    a = Set(split(a, " "))
    b = Set(split(b, " "))
    sort([i for i in union(a, b)]), sort([i for i in intersect(a, b)])
end

sentence1 = "We are really pleased to meet you in our city"
sentence2 = "The city was hit by a really heavy storm"
u, i = func9(sentence1, sentence2)
println(u)
println(i)

SubString{String}["The", "We", "a", "are", "by", "city", "heavy", "hit", "in", "meet", "our", "pleased", "really", "storm", "to", "was", "you"]
SubString{String}["city", "really"]

10. Prime Numbers Array

Given k numbers which are less than n, return the set of prime number
among them Note: The task is to write a program to print all Prime
numbers in an Interval. Definition: A prime number is a natural number
greater than 1 that has no positive divisors other than 1 and itself.

参照先の解答例は n 以下の素数リストを求めるだけのものだけど,ちょっと違うんじゃないか?

function func10(k)
    n = maximum(k)
    primes = [2]
    for i = 3:2:n
        isprime = true
        for j = 2:floor(Int, sqrt(i))
            if mod(i, j) 0
                isprime = false
                break
            end
        end
        if isprime
            append!(primes, i)
        end
    end
    sort([s for s in intersect(Set(k), Set(primes))])
end

n 以下の値の k 個の数のリストから,素数のみから成るリストを求めよ

k = [2, 1, 3, 4, 3, 7, 13, 18, 27, 35]
func10(k)

  2
  3
  7
 13

コメント

Julia で R の deparse() みたいなもの(その2)

2021年01月23日 | ブログラミング

やはり,前の記事のは無理があるので,すっきりさせようと。

ただ,書式を与えられるようにと。

普通にやると動かないのは @printf() がマクロとはいえ,コンパイルされる関数で,fmt の部分に文字変数を指定しても動かない。それを回避するには @eval() を使うというヒントのページの助けを借りて以下のようなものを作った。

julia> using Printf

julia> function dumplist(name::String, object::Array{Float64,2}; fmt="%10.5f")
           nr, nc = size(object)
           println("----- $name -----")
           for i = 1:nr
               for j = 1:nc
                   @eval(@printf($fmt, $(object[i, j])))
               end
               print("\n")
           end
       end
dumplist (generic function with 1 method)

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

julia> dumplist("name", x, fmt="%6.3f")

----- name -----
 1.000 2.000 3.000
 4.000 5.000 6.000

julia> dumplist("name", x, fmt="%6.3g")
----- name -----
     1     2     3
     4     5     6

出力が,パラパラパラ...と出てくるのが可笑しい。

コメント

Julia で R の deparse() みたいなもの

2021年01月22日 | ブログラミング

ちょっと探しても見つからないし,Julia 自体が,変数のタイプがなんなのかは typeof() でわかるんだけど,それに応じて別々の処理をしようとしてもできない。

if typeof(x) == String なんてのができないということ。

そこで,こねくり回して以下のような関数をでっち上げた。

主たる目的は,行列の場合には書式に従って,変数名も一緒にして,綺麗に出力するということだ。変数名を文字列で渡す。

using Printf

function dumplist(argument::String)
    ndims(arg::String) = true
    value = eval(Meta.parse(argument))
    try
        nr, nc = size(value)
    catch
        println("$argument = $value")
        return
    end
    println("----- $argument -----")
    nr, nc = size(value)
    for i = 1:nr
        for j = 1:nc
            @printf("%10.5f", value[i, j])
        end
        print("\n")
    end
end

julia> a = [1 2;3 4]; dumplist("a")
----- a -----
   1.00000   2.00000
   3.00000   4.00000

julia> b = 123.2; dumplist("b")
b = 123.2

julia> c = "string"; dumplist("c")
c = string

julia> d = [1,2,3]; dumplist("d")
d = [1, 2, 3]

 

コメント

Julia で順列,組み合わせなど

2021年01月19日 | ブログラミング

バイナリ・コンビネーション

R の e1071::bincombinatons() に相当するもの。

絶対にお勧めではないが,単純で忘れにくいもの。小さな場合に。

julia> [[x,y,z] for x=0:1 for y=0:1 for z=0:1]
 [0, 0, 0]
 [0, 0, 1]
 [0, 1, 0]
 [0, 1, 1]
 [1, 0, 0]
 [1, 0, 1]
 [1, 1, 0]
 [1, 1, 1]

多分,実際にメモリ内に確保してしまうもの。

julia> all_perm(x, n) = vec(map(collect, Iterators.product(ntuple(_ -> x, n)...)))

julia> all_perm([0, 1], 4)
 [0, 0, 0, 0]
 [1, 0, 0, 0]
 [0, 1, 0, 0]
 [1, 1, 0, 0]
   :略
 [0, 1, 1, 1]
 [1, 1, 1, 1]

julia> for i in all_perm(["a", "b"], 3)
           println(i)
       end
["a", "a", "a"]
["b", "a", "a"]
["a", "b", "a"]
["b", "b", "a"]
["a", "a", "b"]
["b", "a", "b"]
["a", "b", "b"]
["b", "b", "b"]

たぶん,iterator を使って上手くやるもの。

julia> all_perm(x, n) = Iterators.product([x for i = 1:n]...)

julia> for i in all_perm([0, 1], 3)
           println(i)
       end
(0, 0, 0)
(1, 0, 0)
(0, 1, 0)
(1, 1, 0)
(0, 0, 1)
(1, 0, 1)
(0, 1, 1)
(1, 1, 1)

組み合わせ

R の combn() に相当するもの。

julia> using Combinatorics

julia> for i in combinations(1:4, 2)
           println(i)
       end
[1, 2]
[1, 3]
[1, 4]
[2, 3]
[2, 4]
[3, 4]

julia> for i in combinations(["one", "two", "three", "four"], 2)
           println(i)
       end
["one", "two"]
["one", "three"]
["one", "four"]
["two", "three"]
["two", "four"]
["three", "four"]

順列

R の permutations() に相当するもの。

多分,メモリに確保してしまうもの。

julia> for i in permutations([1,2,3])
           println(i)
       end
[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]

多分,iterator で上手くやってくれるもの。

julia> a = [1,2,3]

julia> for i in 1:factorial(length(a))
           println(nthperm(a, i))
       end
[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]

nthperm!() は破壊的なので,全ての順列を網羅できない。
網羅しなくても良い場合や,重複しても良い場合はこちらが多分よいはず。

julia> a = [1,2,3]

julia> for i in 1:factorial(length(a))
           println(nthperm!(a, i))
       end
[1, 2, 3]
[1, 3, 2]
[3, 1, 2]
[1, 2, 3]
[3, 1, 2]
[2, 1, 3]

コメント