RSS

月別アーカイブ: 12月 2021

Julia の勉強:ブラウザに計算結果を表示 (Mux)

Julia 単体ではコンソールアプリケーションしか作れませんが、サードパーティのパッケージを使うことでウェブアプリケーションを作ることができます。

そのひとつに WebIO + Mux を使う方法があります。

そして、WebIO から見て Muxがフロントエンドになります。

これらのパッケージを使うには Julia の対話モードから “]” を押してパッケージモードにし、add コマンドでインストールします。

pkg> add WebIO
pkg> add Mux

次のサンプルは ”Hello, World!” です。(Mux ページのサンプルを参考にしたもの)

using Mux

@app test = (
  Mux.defaults,
  page(respond("<h1>Hello World!</h1>")),
  page("/user/:user", req -> "<h1>Hello, $(req[:params][:user])!</h1>"),
  Mux.notfound()
)

serve(test, 2333)
println("Start 'test' server > ")
read(stdin, Char)

これを実行してブラウザで開くと下のようになります。(アドレスバーに http://localhost:2333 を入力する)

ブラウザのアドレスバーに http://localhost:2333/user/Julia と入力すると、下のように表示されます。

このプログラムの最後で read(stdin, Char) を実行していますが、これがないとプログラムが直ちに終了してしまうので注意してください。

もう少し、複雑な処理をしようとすると Mux 単体ではちょっとたいへんです。

というのも、Mux は HTTP の薄いラッパーで機能があまり高くありません。

それを補うために HttpCommon をインストールしておく必要があります。このパッケージはウェブでよくある共通処理を実装しています。

フォームの例

次の例は、フォームのテンプレートを読み込んでクライアントへ送り、クライアントから POST メソッドで送り返されたフォーム内容を解析して結果をフォームに埋め込んで返すものです。

この例は、2つのベクトルを受け取り、内積を計算し結果をフォームに埋め込んでクライアント (ブラウザ) に返すものです。

起動コマンド
$ julia calc.jl
Start 'calc' server >
GET /
POST /
.
$

Julia ソース (calc.jl)

using Mux
using HttpCommon
using LinearAlgebra

calc_html = "<h1>Error</h1>"

function calc_action(req)
    open("./calc.html", "r") do f
        global calc_html = read(f, String)
    end
    if req[:method] == "POST"
        println("POST /")
        data = req[:data]
        n = length(data)
        param = ""
        for i = 1:n
            param = string(param, Char(data[i]))
        end
        params = parsequerystring(param)  # HttpCommon
        sv1 = split(params["v1"], ",+")  # 空白が + に置き換わっているため
        sv2 = split(params["v2"], ",+")
        v1 = Vector{Int32}(undef, length(sv1))
        v2 = Vector{Int32}(undef, length(sv2))
        n = length(v1)
        for i = 1:n
            v1[i] = parse(Int32, sv1[i])
        end
        n = length(v2)
        for i = 1:n
            v2[i] = parse(Int32, sv2[i])
        end
        # println(v1, "    ", v2)
        y = dot(v1, v2)
        replace(calc_html, "%result%" => y)
    else
        println("GET /")
        replace(calc_html, "%result%" => "")
    end
end

@app calc = (
    Mux.defaults,
    page("/", req -> calc_action(req)),
    Mux.notfound()
)

serve(calc, 2333) # http://localhost:2333
println("Start 'calc' server > ")
read(stdin, Char)

テンプレートファイル (calc.html)

<!DOCTYPE html>
<head>
    <meta name="viewport" content="width=device-width,initial-scale=1">
    <meta charset = "utf-8" />
    <title>Mux test</title>
    <style>
        body {
            margin-left:5%;
            margin-right:5%;
        }

        h1 {
            text-align:center;
            color:crimson;
            padding:8px;
        }

        input[type="text"] {
            border-radius: 3px;
            height:32px;
            padding:2px;
        }

        input[type="button"] {
            border-radius: 3px;
            height:32px;
        }

        fieldset {
            border-width:0px;
        }

        .result {
            font-size:12pt;
            color:magenta;
            padding:5px;
        }
    </style>
    <script>
        function button1Click() {
            if (form1.v1.value == "" || form1.v2.value == "") {
                alert("ベクトルの値が空欄です。");
                return;
            }
            form1.submit();
        }
    </script>
</head>

<body>
    <h1>ベクトルの内積</h1>
    <br />
    <form id="form1" method="POST" action="/" style="margin-left:10%;">
        <p>数値の間は必ず ", " (カンマと空白1つ) にしてください。</p>
        <fieldset>v1: <input type="text" name="v1" id="v1" size="50" placeholder="3, -1, 6" /></fieldset>
        <fieldset>v2: <input type="text" name="v2" id="v2" size="50" placeholder="2, 1, 4" /></fieldset>
        <fieldset><input type="button" id="button1" value=" 計算する " onclick="button1Click()" /></fieldset>
        <br />
        <p class="result" id="result">結果: %result%</p>
    </form>

</body>
</html>

行列式の計算

この例はテンプレートやフォームを使わず、計算結果をブラウザに表示するだけのものです。

具体的には 3×3 の正方行列を乱数 (0:9) を使って10個作り、その行列式を計算して表示しています。

コンソールアプリケーションでは、計算結果が多いと見にくいのでこういう使い方が向いています。

Julia ソース (det.jl)

using LinearAlgebra
using Mux
using Printf

# 正方行列の次元
const N = 3

# 計算とHTML作成
function calc_dets()
    s = "<h1>行列式の計算</h1>\n"
    s *= "<hr /><br />\n"
    for i = 1:10
        A = rand(Int32(0):Int32(9), (N, N))
        d = det(A)
        s *= "<table>"
        for j = 1:N
            s *= "<tr>"
            for k = 1:N
                s *= "<td>"
                if k == 1
                    s *= "| "
                end
                s *= @sprintf("%d", (A[j, k]))
                if k == N
                    s *= " |"
                end
                s *= "</td>"
            end
            s *= "</tr>"
        end
        s *= "</table>\n"
        s *= "<p>= $d</p><br />\n"
    end
    return s
end

# HTML ヘッダ部分
head = """
  <html>
  <head>
  <meta name="viewport" content="width=device-width,initial-scale=1">
  <meta charset = "utf-8" />
  <title>Determinants</title>
  <style>
  body { padding:50px; margin-left:5%; margin-right:5%;}
  h1 { padding:10px; color:crimson; text-align:center;}
  table, td { border-width: 0px; padding:3px;}
  </style>
  </head>
  <body>
"""
# HTML フッタ部分
foot = """
  <p>&nbsp;</p><p>&nbsp;</p>
  </body>
  </html>
"""

# HTML の内容
content = head * calc_dets() * foot

# アプリケーションの定義
@app the_app = (
  Mux.defaults,
  page(respond(content)),
  Mux.notfound()
)

# 実行開始
serve(the_app, 2333)
println("Start 'det' server > ")
read(stdin, Char)
 
コメントする

投稿者: : 2021/12/14 投稿先 Julia

 

タグ: , , , ,

Julia の勉強:平均や標準偏差

Statistics パッケージをインポートすると、平均や標準偏差を簡単に計算できます。

平均

平均は mean() 関数で計算します。標本に欠損値 (missing) があると結果は missing になりますが、skipmissing() 関数を使うと欠損値を避けて平均を計算します。

julia> using Statistics
julia> mean(1:20)
10.5
julia> mean([1, missing, 3])
missing
julia> mean(skipmissing([1, missing, 3]))
2.0

標準偏差

標準偏差は std() 関数で計算できます。

julia> std(1:10)
3.0276503540974917
julia> std([0.5, 0.1, 0.9, 0.3])
0.34156502553198664

分散

分散は var() 関数で計算できます。

julia> var(1:10)
9.166666666666666
julia> var([0.5, 0.1, 0.9, 0.3])
0.11666666666666668

以上は基本的な平均、標準偏差、分散を計算する関数ですが、他にもいくつかのバリエーションがあります。

 
コメントする

投稿者: : 2021/12/11 投稿先 Julia

 

タグ: , , ,

Julia の勉強:乱数

一般的な乱数

0と1 の間の値を取る乱数は Base.rand() で取得できます。

この関数はデフォルトでは上記のような乱数を返しますがオプションを指定すると乱数の範囲、型(整数など)の指定、次元(行列など)の指定もできます。

julia> rand(1:100)  # 整数 1~100の乱数
29
julia> rand(1:100)
34
julia> rand('A':'Z')  # 文字の乱数
'R': ASCII/Unicode U+0052 (category Lu: Letter, uppercase)
julia> rand('A':'Z')
'A': ASCII/Unicode U+0041 (category Lu: Letter, uppercase)
julia> rand(3, 3)  # 3x3の0~1のFloat64
3×3 Matrix{Float64}:
 0.56848   0.138164  0.0733432
 0.539703  0.83782   0.317063
 0.036541  0.623098  0.390291
ulia> rand(Int32, 3, 3)  # Int32の3x3行列
3×3 Matrix{Int32}:
 -1057114174   -999388178   1296397317
 -2070187498   -121350132  -1971020113
   627422604  -1030086586   -361264995

Base.rand() の定義は次のようになっています。

rand([rng=GLOBAL_RNG], [S], [dims...])

rng はシード値、S は型やレンジ、dims は次元です。シード値は次のようにして取得できます。

なお、rand() は Base パッケージに含まれているので、パッケージインポートは不要ですが、それ以外の関数は Random パッケージの関数なので “Random” のインポートが必要です。

julia> using Random
julia> rng = MersenneTwister(1234)

rand() は Base パッケージに含まれるのでインポートしなくても使えます。Random パッケージをインポートすると rand() 以外にも多くの乱数発生機能を利用できます。

次の rand!() はパッケージ Random に含まれる関数です。これは既存の配列を乱数で埋めます。

julia> A = Matrix{Int8}(undef, 2, 2)
2×2 Matrix{Int8}:
 80   34
 32  -82
julia> rand!(A)
2×2 Matrix{Int8}:
 84  -109
 -6   121
julia> rand!(A)
2×2 Matrix{Int8}:
 83  -30
  5   15

正規分布や指数分布を持つ乱数

randn() や randn!() は平均が 0 で標準偏差が 1 の正規分布に従う乱数を発生します。また、randexp() や randexp!() は尺度 1 の指数分布に従う乱数を発生します。

julia> rand()
0.6056531121505694
julia> rand()
0.6306483941786045
julia> rand(1.0:10.0)
6.0
julia> rand(1.0:10.0)
4.0
julia> randexp()
0.1083628567822738
julia> randexp(3)
3-element Vector{Float64}:
 0.11316847843214593
 0.45746821083363404
 0.562927391348907

シャッフル

shuffle() 関数は配列の要素をランダムに入れ替えます。

julia> a = [1, 2, 3, 4, 5]
5-element Vector{Int64}:
 1
 2
 3
 4
 5
julia> shuffle(a)
5-element Vector{Int64}:
 3
 5
 4
 1
 2

部分列や置換

部分列や置換を行う関数がいくつか用意されています。一例としてrandsubseq() を挙げます。

この例では 1:10 のうちから5割の確率で値を取り出してベクトルとして返します。

julia> randsubseq(1:10, 0.5)
7-element Vector{Int64}:
  3
  5
  6
  7
  8
  9
 10

次のように値の組をレンジでなくベクトルでも指定できます。

julia> randsubseq([1, 2, 3, 4, 5], 0.5)
2-element Vector{Int64}:
 1
 3

確率を 1.0 にすると常にすべての要素が選ばれます。

julia> randsubseq([1, 2, 3, 4, 5], 1.0)
5-element Vector{Int64}:
 1
 2
 3
 4
 5

 
コメントする

投稿者: : 2021/12/11 投稿先 Julia

 

タグ: , ,

Julia の勉強: C 言語の関数を使う

Julia では C や Fortran の関数を簡単に利用できます。

Fortran はなじみが薄い言語かもしれませんが、昔は科学技術分野でよく使われていたためライブラリが豊富です。

ここでは、なじみのある C 言語の関数を Julia から呼び出して使ってみます。

動作環境は Linux とします。Windows でも動くかもしれませんが、C ライブラリの形式が違ったりするので Julia のドキュメントのサンプルはそのまま動く可能性は低いです。

C の関数を呼び出すには ccall() を使います。これは関数のように見えますが、「キーワード」になっています。そして、形式としては次の3パターンがあります。

ccall((function_name, library), returntype, (argtype1, ...), argvalue1, ...)
ccall(function_name, returntype, (argtype1, ...), argvalue1, ...)
ccall(function_pointer, returntype, (argtype1, ...), argvalue1, ...)

C 標準でない独自ライブラリ等を使う場合は最初の形式、C 標準ライブラリを使う場合は2番目の形式を使います。3番目は関数ポインタを通して関数を呼び出す場合に使います。

一番簡単な例で clock() (プログラム起動からの経過時間) を呼び出すと次のようになります。

julia> ccall(:clock, Int32, ())
757151

次に初等数学関数 exp(x) を呼び出してみます。おっと、エラーが出てしまいました。これは3番目のパラメータが (Float64) になっているためです。このパラメータはタプルである必要がありますが、(Float64) == Float64 です。

タプルにするには (Float64,) のように “,” が必要です。

julia> ccall(:exp, Float64, (Float64), 1.0)
ERROR: syntax: ccall argument types must be a tuple; try "(T,)" around REPL[2]:1
Stacktrace:
 [1] top-level scope
   @ REPL[2]:1

julia> ccall(:exp, Float64, (Float64,), 1.0)
2.718281828459045

次に独自の C ライブラリを作って Julia から呼び出してみます。このプログラムのファイル名を “addmult.c” とします。

double add(double x, double y) {
  return x + y;
}

double mult(double x, double y) {
  return x * y;
}

double inv(double x) {
  return -x;
}

これをビルドして .so ファイルにします。

$ gcc -shared -fPIC -o addmult.so addmult.c

念のために C のテストプログラムを作ってちゃんと動くことを確認します。

#include <stdio.h>

double add(double, double);
double mult(double, double);
double inv(double);

int main() {
  double y;
  y = add(1.5, 5.5);
  printf("%f\n", y);
  y = mult(2.0, 4.0);
  printf("%f\n", y);
  y = inv(1.0);
  printf("%f\n", y);
  return 0;
}

これをビルドして動作確認を行います。

$ gcc test_addmult.c -o test_addmult ./addmult.so
$ ./test_addmult
7.000000
8.000000
-1.000000

いよいよ Julia からこの C ライブラリの関数を呼び出してみます。

y = ccall((:add, "./addmult.so"), Float64, (Float64, Float64), 1.5, 5.5)
println(y)
y = ccall((:mult, "./addmult.so"), Float64, (Float64, Float64), 2.0, 4.0)
println(y)
y = ccall((:inv, "./addmult.so"), Float64, (Float64,), 1.0)
println(y)

このプログラムの名前は call_addmult.jl です。下のようにこのプログラムを julia で実行します。

$ julia call_addmult.jl
7.0
8.0
-1.0

期待通りに結果が表示されましたね。

 
コメントする

投稿者: : 2021/12/10 投稿先 Julia

 

タグ: , ,

Julia の勉強:スレッド

Julia では並列処理のためのメカニズムとしてタスク以外にスレッドと分散計算をサポートしています。

タスクは関数ごとに設定可能でしたが、スレッドはもっと細かくてある処理の中で並列処理を可能にします。

Julia の対話モードでマルチスレッドを扱うには起動時に -t スイッチでスレッド数を指定する必要があります。(Julia 1.4 以前では環境変数 JULIA_NUM_THREADS にスレッド数を設定しておく必要があります)

スレッド数は Threads.nthreads() 関数で確認できます。

user@raspi4:~$ julia -t 4
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.0 (2021-11-30)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> Threads.nthreads()
4

スレッドは @threads マクロを使うと簡単に利用できます。このマクロは並列ループをサポートします。つまり、このマクロで作ったループは並列に実行されます。

# スレッドテスト "julia -t 4" で実行すること。
a = Vector{Int64}(undef, 10)
Threads.@threads for i = 1:10
    a[i] = Threads.threadid()
end
println(a)

実行例

 julia -t 4 thread_test.jl
[1, 1, 1, 2, 2, 2, 3, 3, 4, 4]

このサンプルではベクトル a の各要素にスレッド ID が並列ループ内で書き込まれ10回繰り返すトループを抜けてベクトルの内容が表示されます。

 
コメントする

投稿者: : 2021/12/09 投稿先 Julia

 

Julia の勉強:タスク

タスクは @task マクロで作ります。これは Task(()->x) と同じです。タスクの実行はタスクを作っただけでは始まりません。schedule() 関数を実行する必要があります。

julia> t = @task begin; sleep(5); println("done"); end
Task (runnable) @0x000000000ed7a0e0
julia> schedule(t)
Task (runnable) @0x000000000ed7a0e0
julia> done
julia>
julia> t = Task(()->begin; sleep(5); println("done"); end)
Task (runnable) @0x000000000f3d9140
julia> schedule(t)
Task (runnable) @0x000000000f3d9140
julia> done
julia>

wait() を使うと他のタスクが終了するまでブロックできます。

次の例では、タスク t が終了すると他のタスクがないため直ちに “julia>” が表示されます。

julia>  t = @task begin; sleep(5); println("done"); end
Task (runnable) @0x000000000f4c6720
julia> schedule(t); wait(t)
done
julia>

次の例ではタスク t, u を同時に実行し、時間がかかる方のタスク u が終わってからプログラムを終了します。

# 2つのタスク
t = @task begin; sleep(5); println("done t"); end
u = @task begin; sleep(6); println("done u"); end
schedule(t);schedule(u);wait(u)  # タスク u が終わるまで待つ。

タスクの生成と開始は同時に使われるので、便利なマクロ @async が用意されています。@async は schedule(@task x) と同じです。

# 2つのタスク (@sync)
t = @async begin; sleep(5); println("done t"); end
u = @async begin; sleep(6); println("done u"); end
wait(u)  # タスク u が終わるまで待つ。

タスク t が終了したかどうかは istaskdone(t) 関数で取得できます。

wait(t) の代わりに fetch(t) を使うと、タスクで実行した関数の戻り値を取得できます。

# 関数の戻り値
function func(n)
    nn = 1
    for i = 1:n
        nn *= i
    end
    return nn
end
t = @async(func(4))
nn = fetch(t)
println(nn)

このプログラムを実行すると 24 と表示されます。

 
コメントする

投稿者: : 2021/12/08 投稿先 Julia

 

タグ: ,

Julia の勉強:TOML 形式のファイルを読む

Julia では設定ファイルとして TOML 形式のファイルがよく使用されるため、TOML パッケージが標準として含まれています。

TOML ファイルは “INI” ファイルによく似ています。

そのため、ちょっとした設定ファイルとしては見やすく簡潔です。

基本的には TOML ファイルを読んでその内容を parse() で解析すれば辞書としてその内容を得ることができます。

# config.toml を読む。
using TOML
open("config.toml", "r") do io
    s = read(io, String)
    #println(s)
    conf = TOML.parse(s)
    println("Basic.title: ", conf["Basic"]["title"])
    println("Basic.date: ", conf["Basic"]["date"])
    println("Basic.version: ", conf["Basic"]["version"])
    println("Description.content: ", conf["Description"]["content"])
end

ここで “config.toml” は下のようになっています。

[Basic]
  title = "TOML Test"
  date = "2021-12-01"
  version = "1.0.1"
[Description]
  content = "TOML のテストです。"

このプログラムは parsefile() 関数を使えばもっと簡単になります。

# config.toml を読む。
using TOML
conf = TOML.parsefile("config.toml")
println("Basic.title: ", conf["Basic"]["title"])
println("Basic.date: ", conf["Basic"]["date"])
println("Basic.version: ", conf["Basic"]["version"])
println("Description.content: ", conf["Description"]["content"])

このプログラムを実行すると下のようになります。

Basic.title: TOML Test
Basic.date: 2021-12-01
Basic.version: 1.0.1
Description.content: TOML のテストです。

 
コメントする

投稿者: : 2021/12/08 投稿先 Julia

 

タグ: ,

Julia の勉強:線形代数 (LinearAlgebra)

Julia には LinearAlgebra という「標準パッケージ」があります。

標準とはいっても Core や Base パッケージのように組み込みではないのでインストールが必要です。Julia の対話モードから “]” を押して、パッケージモードへ移行し次のコマンドを実行します。

pkg> add LinearAlgebra

線形代数の勉強お勧め動画 (予備校のノリで学ぶ線形代数入門)

  1. 概観&ベクトル
  2. 行列 (演算規則)
  3. 一次変換
  4. 一次独立と一次従属
  5. 連立方程式 (掃き出し法)
  6. 不定と不能
  7. 階数 (ランク)
  8. 行列式 (定義と性質)
  9. 余因子展開
  10. 逆行列:定義
  11. 逆行列:掃き出し法
  12. 固有値・固有ベクトル
  13. 対角化 (重解がない場合)
  14. 対角化 (重解がある場合)

ベクトル

ベクトルの初期化

横ベクトル は v1 = [0.5, 1.1, -2.0] 、縦ベクトルは v2 = [3.3; 1.7; 0.2] のようにします。要素の型を指定する場合は次のようにします。

julia> v2 = Int32[1, 2, 4]
3-element Array{Int32,1}:
 1
 2
 4

要素の値が未定義のベクトルは v3 = Vector{Float64}(undef, 3) のようにします。undef の代わりに nothing や missing も使えます。

julia> v3 = Vector{Float64}(undef, 3)
3-element Array{Float64,1}:
 1.390662761733164e-309
 1.390662959311913e-309
 1.39066270598319e-309

ベクトルの絶対値 (ノルム)

ベクトルの絶対値ですが、abs(v) とやるとエラーになります。norm(v) が正しいです。

julia> v2
2-element Vector{Int64}:
 1
 2
julia> norm(v2)
2.23606797749979
julia> sqrt(5.0)  # 1.0^2 + 2.0^2 = 5.0
2.23606797749979

内積と外積

次にベクトルの内積ですが、v1 * v2 とやるとエラーになります。dot(v1, v2) が正しいです。また、外積は cross(v1, v2) で計算できます。

using LinearAlgebra
v1 = [0.5, 1.1, -2.0]
v2 = [3.3; 1.7; 0.2]
# 絶対値
println(norm(v1))
# 内積
println(dot(v1, v2))
# 外積
println(cross(v1, v2))

実行結果

2.3366642891095846
3.12
[3.62, -6.699999999999999, -2.78]

行列

初期化

行列を初期化する場合、行と行の間の区切りに “;” を使います。同じ行の要素の区切りは ” ” を使います。

julia> A = [3 1 3; 0 -1 1; 5 2 -2]
3×3 Array{Int64,2}:
 3   1   3
 0  -1   1
 5   2  -2

行列の型を特に指定しないと64ビット整数または64ビット浮動小数点数になりますが、型を指定したい場合は次のよう行列の前に型名を指定します。

julia> A = Int32[3 1 3; 0 -1 1; 5 2 -2]
3×3 Array{Int32,2}:
 3   1   3
 0  -1   1
 5   2  -2

要素を後から決める場合は、次のように要素が undef の行列を作ります。

julia> A = Matrix{Float64}(undef, 3, 3)
3×3 Array{Float64,2}:
 5.0e-323      6.0e-323      6.4e-323
 1.0e-323      1.0e-323      1.0e-323
 1.39066e-309  1.39066e-309  1.39066e-309

要素に値を代入するには次のように行います。ここで先頭の要素のインデックスは [0, 0]ではなく [1, 1]

になります。

julia> A[1, 1] = 0.0
0.0
julia> A
3×3 Array{Float64,2}:
 0.0           6.0e-323      6.4e-323
 1.0e-323      1.0e-323      1.0e-323
 1.39066e-309  1.39066e-309  1.39066e-309

インデックスの指定ですが、最初の方が行番号、後の方が列番号になります。

julia> A[2, 1] = 5.0
5.0
julia> A
3×3 Array{Float64,2}:
 0.0           6.0e-323      6.4e-323
 5.0           1.0e-323      1.0e-323
 1.39066e-309  1.39066e-309  1.39066e-309

単位行列を作るには次のようにします。

julia> E = Matrix{Int32}(I, 3, 3)
3×3 Array{Int32,2}:
 1  0  0
 0  1  0
 0  0  1

要素がすべて 0 の行列は次のように作ります。

julia> Z = zeros(Int32, 3, 3)
3×3 Array{Int32,2}:
 0  0  0
 0  0  0
 0  0  0

要素がすべて 1 の行列も次のようにして作れます。

julia> O = ones(Int32, 3, 3)
3×3 Array{Int32,2}:
 1  1  1
 1  1  1
 1  1  1

乱数で初期化することもできます。下の例ではシード値等はデフォルト値を使っていますが、もっと細かい指定もできます。

julia> R =rand(3, 3)
3×3 Array{Float64,2}:
 0.588734  0.609012  0.680494
 0.212418  0.1448    0.87533
 0.380341  0.219114  0.570458

行列の演算

行列の演算は足し算、引き算、掛け算がありますが、割り算はありません。ベクトルの掛け算 (内積) では “*” は使えませんでしたが、行列では “*” が使えます。

また、ベクトルも行列の一種なので、行列とベクトルの掛け算も可能です。

足し算、引き算は同じ次元の行列同士でのみ可能です。掛け算は最初の行列の列数と後の行列の行数が同じである必要があります。

次のサンプルは行列の簡単な演算例です。

# 行列の演算
using LinearAlgebra
# 定義
A = Int32[5 9 1; 3 -1 -2; 1 0 8]
println("A = ", A)
B = Int32[1 5 -1; -8 4 4; 2 4 -8]
println("B = ", B)
C = Int32[-5 5 -9; 1 -6 6; -1 1 7]
println("C = ", C)
E = Matrix{Int32}(I, 3, 3)
println("E = ", E)
x = Int32[2; 1; -6]
println("x = ", x)
# 計算
println("A + B = ", A + B)
println("A - B = ", A - B)
println("A * B = ", A * B)
println("A * x = ", A * x)
println("+ 結合法則 ", (A + B) + C == A + (B + C))
println("* 結合法則 ", (A * B) * C == A * (B * C))
println("+ 交換法則 ", A + B == B + A)
println("* 交換法則 ", A * B == B * A)  # 等しくない
println("* 分配法則 (スカラー) ", 2.0 * (A + B) == 2.0 * A + 2.0 * B)
println("* 分配法則 (行列) ", C * (A + B) == C * A + C * B)
println("* 単位行列 (交換法則) ", E * A == A * E)
# 正方行列以外の掛け算
D = Int32[2 1 0; 0 1 2]
println("D = ", D)
F = Int32[1 2; 3 4; 5 6]
println("F = ", F)
println("D * F = ", D * F)
# べき乗
println("A^2 = ", A^2)
# スカラー倍
println("2 * A = ", 2 * A)0 * A)

実行結果

A = Int32[5 9 1; 3 -1 -2; 1 0 8]
B = Int32[1 5 -1; -8 4 4; 2 4 -8]
C = Int32[-5 5 -9; 1 -6 6; -1 1 7]
E = Int32[1 0 0; 0 1 0; 0 0 1]
x = Int32[2, 1, -6]
A + B = Int32[6 14 0; -5 3 2; 3 4 0]
A - B = Int32[4 4 2; 11 -5 -6; -1 -4 16]
A * B = Int32[-65 65 23; 7 3 9; 17 37 -65]
A * x = Int32[13, 17, -46]
+ 結合法則 true
* 結合法則 true
+ 交換法則 true
* 交換法則 false
* 分配法則 (スカラー) true
* 分配法則 (行列) true
* 単位行列 (交換法則) true
D = Int32[2 1 0; 0 1 2]
F = Int32[1 2; 3 4; 5 6]
D * F = Int32[5 8; 13 16]
A^2 = Int32[53 36 -5; 10 28 -11; 13 9 65]
2 * A = [10 18 2; 6 -2 -4; 2 0 16]

一次変換

行列にベクトルを掛けるということは、そのベクトルを別のベクトルに変換するということです。この変換のことを一次変換と言います。(一次式だけからなるため)

単位行列にベクトルを掛けてもは元のベクトルのままです。

julia> E = Matrix{Float32}(I, 3, 3)
3×3 Matrix{Float32}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
julia> x = [2.0; 5.0; 1.0]
3-element Vector{Float64}:
 2.0
 5.0
 1.0
julia> E * x
3-element Vector{Float64}:
 2.0
 5.0
 1.0

対角行列にベクトルを掛けると、ベクトルの方向は変わらず長さだけが変わります。

julia> D = Float32[2.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 3.0]
3×3 Matrix{Float32}:
 2.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  3.0
julia> D * x
3-element Vector{Float64}:
 4.0
 5.0
 3.0

次のように定義される回転行列を掛けるとベクトルの偏角が変わりベクトルを回転させることができます。

julia> theta = pi / 4.0
0.7853981633974483
julia> T = [cos(theta) -sin(theta); sin(theta) cos(theta)]
2×2 Matrix{Float64}:
 0.707107  -0.707107
 0.707107   0.707107
julia> v = [1.0; 1.0]
2-element Vector{Float64}:
 1.0
 1.0
julia> T * v
2-element Vector{Float64}:
 1.1102230246251565e-16
 1.414213562373095

対角行列と回転行列を組み合わせてベクトルに掛けると長さと角度を変えたベクトルに変換できます。

これらは二次元だけでなく n 次元で定義されています。

一次独立と一次従属

n 個のベクトルの集合 (v1, v2, … vn) の中の任意のベクトル vk が残りすべてのベクトルの線形結合で表すことができなければ「一次独立」と言い、そうでなければ「一次従属」と言います。

線形結合とは c1*v1 + c2*v2 + … + cn*vn のように定数 ck とベクトル vk を掛けたものの和です。

2次元で言えば、一次独立は方向の異なるベクトルどうし、一次従属は方向が同じなベクトルどうしになります・

つまり、v2 = c1 * v1 のような関係です。3次元で言えば、3つのベクトルのうち他のベクトルが作る平面上にそのベクトルが存在すれば一次従属です。

つまり、c1*v1 + c2*v2 = v3 になるような数 c1, c2 が存在すれば一次従属です。

一次独立かどうかの判別は「ランク (階数)」によって知ることができます。

つまり、n 個のベクトルの集合 (v1, v2, … vn) のランクが n ならこれらのベクトルは一次独立、n 未満なら一次従属なベクトルが存在することになります。

Julia ではランクは LinearAlgebra.rank() で得ることができます。

次の例で行列 A で定義される3つのベクトルは一次独立ですが、B のランクは3つのベクトルに対してランクが 2 なので一次従属です。

julia> using LinearAlgebra
julia> A = [1 0 0; 0 2 1; 2 1 0]
3×3 Matrix{Int64}:
 1  0  0
 0  2  1
 2  1  0
julia> B = [2 0 0; 1 0 0; 2 1 0]
3×3 Matrix{Int64}:
 2  0  0
 1  0  0
 2  1  0
julia> rank(A)
3
julia> rank(B)
2

連立一次方程式

連立一次方程式を解くには各式の係数と右辺のベクトルを1つにした係数行列を作って「掃き出し法」を使うと機械的に解くことができるのでよく使用されます。

しかし、Julia ではずばり連立方程式を解くための関数があるのでそれを利用します。

ところで、連立方程式は必ず解けるとは限りません。係数行列のすべての行が一次独立でないと解くことができません。

つまり、係数行列のランクが係数行列の次元より小さいと解くことができません。

#  連立方程式
using LinearAlgebra

# Ax = b を LAPACK.gesvx! を使って解く。
A = [0.5 0.4 1.2; -0.8 -03 1.5; 9.1 4.4 -2.0]
b = [0.1; 5.2; -1.0]  # 注意:[0.1 5.2 -1.0] だと Matrix になってしまう。Vector でないとエラーになる。
tx = LAPACK.gesvx!(A, b)
x = tx[1]  # タプルの最初の要素が答え
println(x)
# 検算
y = A * x
println(y)

実行結果

[0.8275653436592451; -1.7869373184898356; 0.33416021297192644]
[0.09999999999999998; 5.2; -0.9999999999999998]

連立方程式の解を求める関数は、行列が特殊な場合などについていろいろ用意されています。

行列式

行列式は正方行列を基にある手順で計算したスカラー値です。

2次と3次の正方行列から計算される行列式の値は簡単に計算できますが、それ以上の正方行列の場合は手計算では面倒です。

Julia ではもちろん行列式の計算も関数を呼び出すことにより簡単に計算できます。

行列式を計算する関数は det() です。

julia> using LinearAlgebra
julia> A = Int32[0 2 2; 1 -1 0; 3 1 2]
3×3 Matrix{Int32}:
 0   2  2
 1  -1  0
 3   1  2
julia> det(A)
4.0

行列式は「余因子展開」で使用されたり、2次元空間での閉じた図形を一次変換したときの面積の倍率に対応しています。

余因子展開

余因子とは

余因子とは行列の m 行 n 列の要素を取り去った小行列の行列式にある決められた方法で符号をつけたものです。

決められた方法とは m + n が偶数なら正、奇数なら負とするものです。

次の例では 3×3 正方行列 A の 3 行 3 列 を除いた 行 1:2, 列 1:2 を取り出した行列 B の余因子行列を計算しています。余因子の符号は 3 + 3 なのでそのままになります。

julia> A = Int32[1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int32}:
 1  2  3
 4  5  6
 7  8  9
julia> B = A[1:2, 1:2]
2×2 Matrix{Int32}:
 1  2
 4  5
julia> det(B)
-3.0

余因子展開

余因子展開は行列式の計算をより低次の行列式に帰着させたものです。3次正方行列で2行2列に着目した場合で言えば、次のようになります。a21 などは行列の要素の値、A21 などは余因子です。

det(A) = a21*​A21​+a22*​A22​+a23*​A23​

逆行列

逆行列は元の行列と掛けると単位行列になるという行列です。Julia では関数 inv() で簡単に逆行列を計算できます。つまり、E を単位行列とすれば次のようになります。

A * inv(A) = E

実際に逆行列を計算してみます。

julia> A = [0.0 2.1 5.2; 3.4 8.2 -1.1; 7.9 -2.5 6.9]
3×3 Matrix{Float64}:
 0.0   2.1   5.2
 3.4   8.2  -1.1
 7.9  -2.5   6.9
julia> B = inv(A)
3×3 Matrix{Float64}:
 -0.120003   0.0612835   0.100207
  0.071672   0.0915797  -0.0394141
  0.163363  -0.0369841   0.0159172
julia> A * B
3×3 Matrix{Float64}:
 1.0           2.77556e-17  -1.38778e-17
 2.77556e-17   1.0           3.46945e-18
 2.22045e-16  -5.55112e-17   1.0

行列の掛け算では交換法則は成り立ちませんが、逆行列との掛け算では交換法則が成り立ちます。

julia> B * A
3×3 Matrix{Float64}:
  1.0          0.0           0.0
  0.0          1.0          -5.55112e-17
 -2.77556e-17  6.93889e-18   1.0

この行列は正確には単位行列ではありませんが、ほぼ単位行列になっています。これはコンピュータによる浮動小数点演算の誤差によるものです。ほぼ等しいかどうかは isapprox() 関数で知ることができます。

julia> A = [1.0 0.0 0.0; 0.0 1.0 -5.55112e-17; -2.77556e-17 6.93889e-18 1.0]
3×3 Matrix{Float64}:
  1.0          0.0           0.0
  0.0          1.0          -5.55112e-17
 -2.77556e-17  6.93889e-18   1.0
julia> E = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
julia> isapprox(A, E)
true

(参考) A*x = B を満たすベクトル x は B \ A で求めることができます。

固有値

固有値とはベクトルを線形変換したとき、元のベクトルのスカラー倍になっているとき、そのスカラー値を固有値と言います。また、そのベクトルを固有ベクトルと言います。

Julia では固有値と固有ベクトルは eigen() 関数で計算できます。

julia> A = [2.0 2.4 1.0; -0.6 -1.7 0.0; 1.8 3.8 2.4]
3×3 Matrix{Float64}:
  2.0   2.4  1.0
 -0.6  -1.7  0.0
  1.8   3.8  2.4
julia> eivals, eivec = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 -1.4137209983476255
  0.8506178385312904
  3.2631031598163345
vectors:
3×3 Matrix{Float64}:
  0.353439  -0.845957   0.540171
 -0.740758   0.199001  -0.0653024
  0.571277   0.494727   0.839018

これの意味ですが、values が固有値 (3つある)、vectors が固有ベクトル (3つある) です。

一つ目の固有ベクトルは [0.353439; -0.740758; 0.571277] です。2つ目、3つ目も同様に3×3 行列の列になります。

実際に元の行列に固有ベクトルを掛けると、固有値 * 固有ベクトルになるか検算してみます。

julia> x = [0.353439; -0.740758; 0.571277]
3-element Vector{Float64}:
  0.353439
 -0.740758
  0.571277
julia> A * x
3-element Vector{Float64}:
 -0.4996642
  1.0472252000000002
 -0.8076253999999996
julia> -1.4137209983476255 * 0.353439
-0.49966413593498643
julia> -0.740758 * -1.4137209983476255
1.0472251392939904
julia> 0.571277 * -1.4137209983476255
-0.8076262907730365

対角化

対角行列とは対角成分以外 0 で、対角成分が正方行列の固有値であるような行列です。

n 次正方行列 A の n 個の一次独立な固有ベクトル x1, x2, … xn を並べた行列を P とすると行列 A は次のようにして対角化できます。

inv(P) * A * P

対角行列ではべき乗が簡単にできます。つまり、対角行列の n 乗は対角成分の n 乗になります。

これを Julia で確認してみます。

julia> A
3×3 Matrix{Float64}:
  2.0   2.4  1.0
 -0.6  -1.7  0.0
  1.8   3.8  2.4
julia> lambda, P = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 -1.4137209983476255
  0.8506178385312904
  3.2631031598163345
vectors:
3×3 Matrix{Float64}:
  0.353439  -0.845957   0.540171
 -0.740758   0.199001  -0.0653024
  0.571277   0.494727   0.839018
julia> D = inv(P) * A * P
3×3 Matrix{Float64}:
 -1.41372       6.31439e-16  -2.22045e-16
 -1.66533e-16   0.850618      2.77556e-16
  0.0          -4.44089e-16   3.2631
julia> D ^ 2
3×3 Matrix{Float64}:
 1.99861      -3.55565e-16  -4.10645e-16
 9.37755e-17   0.723551      1.14179e-15
 7.39557e-32  -1.82686e-15  10.6478

 
コメントする

投稿者: : 2021/12/07 投稿先 Julia

 

タグ: ,

Julia の勉強: Julia 独特の演算子

Julia には他の言語に見られない独特の演算子があってとまどうことがあります。

型の派生

“<:” 演算子は基本の型をベースに型を定義するときに使われます。

abstract type Number end
abstract type Real     <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer  <: Real end
abstract type Signed   <: Integer end
abstract type Unsigned <: Integer end
primitive type Float16 <: AbstractFloat 16 end
primitive type Float32 <: AbstractFloat 32 end
primitive type Float64 <: AbstractFloat 64 end

型の注釈

“::” 演算子は型のヒントをコンパイラに与えるのに使います。Julia では動的に変数の型が決まるのであくまで型をしばるものではなく注釈のようなものです。(Python の型ヒントと同じ)

 struct Foo
           bar
           baz::Int
           qux::Float64
       end

関数の合成

” 演算子は関数のパラメータに何重にも関数を取る場合、それらの関数を外に出して見やすくできるものです。

これは sqrt(3 + 6) と同じ意味です。

julia> (sqrt ∘ +)(3, 6)
3.0

関数のパイプライン

“|>” 演算子はパラメータに対して上流から下流に向かって関数を適用する演算子です。

julia> 1:10 |> sum |> sqrt
7.416198487095663

これは、sqrt(sum(1:10)) と同じです。

ベクトルに全要素に関数を適用する

“.” 演算子を関数パラメータの前に置くと、パラメータがベクトルだった時、そのすべての要素に関数が適用されます。

julia> A = [1.0, 2.0, 3.0]
3-element Vector{Float64}:
 1.0
 2.0
 3.0

julia> sin.(A)
3-element Vector{Float64}:
 0.8414709848078965
 0.9092974268256817
 0.1411200080598672

ベクトル・行列の掛け算、割り算

演算子 “.*”, “./” は要素ごとの掛け算、割り算を行います。

julia> a = [2.0, 4.0]
2-element Vector{Float64}:
2.0
4.0
julia> a ./ 2.0
2-element Vector{Float64}:
1.0
2.0

演算子 “/” は a * inv(b) を計算します。

julia> a = [0.5 0.4; 1.0 -1.0]
2×2 Matrix{Float64}:
 0.5   0.4
 1.0  -1.0

julia> b = [2.0 1.0; 1.5 0.5]
2×2 Matrix{Float64}:
 2.0  1.0
 1.5  0.5

julia> a / b
2×2 Matrix{Float64}:
  0.7  -0.6
 -4.0   6.0

julia> a * inv(b)
2×2 Matrix{Float64}:
  0.7  -0.6
 -4.0   6.0

有理数

“//” を使うと有理数を表せます。例えば、1 // 2 は 2分の1 (= 0.5) を意味します。

julia> 1//2
1//2
julia> 1//2==0.5
true
 
コメントする

投稿者: : 2021/12/04 投稿先 Julia

 

タグ: ,

Julia の勉強:SQLite

Julia で使用するデータをデータベースで管理すると便利なことがあります。

ここでは軽量なデータベースである SQLite3 を Julia で利用してみます。

まず、SQLite パッケージをインストールしておきます。(DataFrames のインストールがまだなら同時にインストールします。)

pkg> add SQLite

このパッケージの詳細は SQLite.jl を参照してください。

以下のサンプルで使用するデータベースは “test.db” というファイルで、使用するテーブルは次のように定義されています。

CREATE TABLE m3x3 (id integer not null primary key autoincrement, row1 text not null, row2 text not null, row3 text not null, info text);

これからわかるようにテーブル “m3x3” は、主キーが “id” で 3×3 行列を定義しています。テーブルのカラムは “row1”, “row2”, “row3”, “info” でいずれも text 型です。

row1, row2, row3 は 3×3 行列の行データで空白区切りで3つの数を格納しています。”info” は備考です。

テーブルを読み込むサンプル

このサンプルはテーブル m3x3 を読んで DataFrame に格納し、その内容を表示します。

表示方法は前半では DataFrame のデータを表示すします。後半では eachrow 関数を使い DataFrame の行を取り出して表示しています。

# Query test.m3x3
using SQLite
using DataFrames

db = SQLite.DB("./test.db")
result = DBInterface.execute(db, "SELECT id, row1, row2, row3 FROM m3x3") |> DataFrame
println(typeof(result))
n = nrow(result)
println("nrow = ", n)
for i = 1:n
  println(result[i, "row1"], ";", result[i, "row2"], ";", result[i, "row3"])
end
# eachrow
println("\neachrow")
rows = eachrow(result)
for r in rows
    println(r.id)
end

テーブルにデータを挿入するサンプル

このサンプルは1行のデータをテーブル “m3x3” に挿入します。

using SQLite
sql = "INSERT INTO m3x3(row1, row2, row3, info) VALUES(3.5, 5.5, 7.5, 'TEST')"
db = SQLite.DB("./test.db")
DBInterface.execute(db, sql)
println("Done.")

 
コメントする

投稿者: : 2021/12/04 投稿先 Julia

 

タグ: ,