kivantium活動日記

プログラムを使っていろいろやります

2018年4月上旬

近況報告

大学院に進学しました。
23歳になりました。

この時期の話題

化合物の合成経路の自動設計

AlphaGoと同じようなMCTSニューラルネットワークを組み合わせる方法で化合物の合成方法を見つけられるという話。

日本人英語 : 英語発音の実態とその矯正法

日本人英語の矯正法について

東大英語部会の発音記号の指導
todai.tv

技術系英語Podcast

技術についてのPodcast。英語の勉強に良いかもしれない?

英語の形容詞の順番

Adjectives: order - English Grammar Today - Cambridge Dictionary

She was a beautiful tall thin young black-haired Scottish woman. というようなたくさんの形容詞が付く場合にちゃんと順番が決まっているよという話。

ターミナルでVimを起動したときの色



テレンス・タオのルベーグ積分本 PDF

公開されているのを知らなかった。

moedict

中国語辞典。特に萌えとは関係ないっぽい。

メンヘラちゃん。(LINEスタンプ)

中国で流行っているらしい。かわいい。

富岡製糸場のブリュナエンジン





その他








読んだ本

百合ドリル (MFC キューンシリーズ)

百合ドリル (MFC キューンシリーズ)

ときめく、はじめての。 (百合姫コミックス)

ときめく、はじめての。 (百合姫コミックス)

2018年3月 第1,2,3週

今週のTwitter


ニュース




スキー旅行



その他

  • ELSAという発音矯正アプリが良さそう

Rosetta Code

プログラミング言語の比較サイト。ロゼッタストーンのように同じ処理を異なる言語で書いて紹介している。
Category:Programming Tasks - Rosetta Codeを見るのが見やすそう。

日本の古本屋

古書店の在庫を横断検索できる。

Blender Guru

Blenderのすごいチュートリアルがいっぱいあるらしい。

量子コンピューティングの講義

読んだ本

三ツ星カラーズ5

最後にして最初のアイドル【短篇版】ひとりぼっちの○○生活(4)
最後にして最初のアイドル【短篇版】

最後にして最初のアイドル【短篇版】

きんいろモザイク アンソロジーコミック (1) 川柳少女(1) われはロボット 〔決定版〕

2018年2月 第3,4週

Twitterより








NSGA-II

多目的最適化という目的関数が複数ある場合の最適化の問題設定があるらしい。単純な線形和で1つの目的関数に落とす方法しか知らなかったので勉強したい。
NSGA-IIという遺伝的アルゴリズムがよく使われるらしい。
qiita.com

Authorea

論文執筆サービス。Overleafより使いやすいという評判を見るが、具体的にどういいのかは試していないので分からない。
www.authorea.com

Hyper Collocation

hypcol.marutank.net

Keywords

読んだ本

私は皆様のおすすめの本・漫画・アニメ・ゲームなどを知りたいです。
自分が好きなコンテンツの話をどんどんして推薦してくれると嬉しいです。

で、私が先週・今週読んだ本の話。

無能なナナ 1-3巻

ネタバレがひどいので、ググったり表紙を見たりせずに公式サイトで公開されている1話をまずは読むべき。

1話の展開が話題になりましたが、2話以降も面白さを維持していてとても良かった。
特に良いのは、(主人公視点の情報は隠すものの)敵視点の情報はきちんと描写した後で事件が起こり、後から整合性のある説明がきちんと与えられるという形式になっていること。
能力バトルだと読者が知り得ない新能力が突然出てくるなどして興ざめすることが多いですが、主人公の性質上そういうことはないため納得できる動機→行動→結果が提示されるのが良い。

極黒のブリュンヒルデ6巻のスカジ(未来予知ができる上に予知の実現のために未来干渉ができる敵)との戦いも、主人公の行動は一部省略され敵目線では完全な描写→事件→伏線を回収する解説の流れが見事でとても好き。

百万光年のちょっと先

冬木糸一さんのレビューを見て買った短編集。
面白いのもいくつかあったのだが、好きでない作品も多かった。


というのがその理由だろう。
現在の科学の常識から想像できるトリック以外を使われると何でも許される気がして興ざめしてしまう。


この発言は、

という感想を念頭に置いて書いたもの。

気に入った短編もあった。



百万光年のちょっと先 (JUMP j BOOKS)

百万光年のちょっと先 (JUMP j BOOKS)

アリスマ王の愛した魔物

小川一水さんは天冥の標シリーズでとても好きになった作家。短編集が出たとのことで買った。
特に太陽系に人類が出たのを前提とした上でどうするかという話が好きで、この短編集でいうと「ゴールデンブレッド」と「星のみなとのオペレーター」が良かった。

ムダヅモ無き改革

全巻5円セールなので買った。麻雀が何も分からない人でも勢いで押し切って読ませるのがすごい。
技術立国日本の描写がASIMO中心なのが今となっては悲しい。

情報幾何学の新展開

Amazonでは品切れなのでサイエンス社のホームページから注文すると良さそう。



動物の賢さがわかるほど人間は賢いのか

「動物の賢さがわかるほど人間は賢くなってきた」という話をするためにたくさん例を上げて説明していた。流し読みしてしまった。

動物の賢さがわかるほど人間は賢いのか

動物の賢さがわかるほど人間は賢いのか

2018年2月第1,2週

Twitter









読んだ本

ゆるゆり: 16 (百合姫コミックス)

ゆるゆり: 16 (百合姫コミックス)

岩波講座 応用数学〈21〉〔対象10〕 システムと制御の数理/〔対象12〕 情報幾何の方法

岩波講座 応用数学〈21〉〔対象10〕 システムと制御の数理/〔対象12〕 情報幾何の方法

RustによるNelder-Mead法の実装

Rustの練習その2としてNelder-Nead法を実装しました。

Nelder-Mead法のアルゴリズムには様々な変種がありますが、ここでは以下の本に書いてあるものを使いました。

Derivative-Free and Blackbox Optimization (Springer Series in Operations Research and Financial Engineering)

Derivative-Free and Blackbox Optimization (Springer Series in Operations Research and Financial Engineering)

また、最適化の対象にはTest functions for optimization - Wikipediaからいくつか選びました。

実装

行列演算ライブラリとしてrust-ndarrayを使うためにCargo.toml

[dependencies]
ndarray = "0.11.0"

を追記しました。

main.rsは以下の通りです。

extern crate ndarray;

use ndarray::*;
use std::f64;
use std::f64::consts;

fn booth(a: &Array1<f64>) -> f64 {
    if a.len() != 2 {
        panic!("input dimension must be 2.");
    }
    let x = a[0];
    let y = a[1];
    (x+2.0*y-7.0).powf(2.0) + (2.0*x+y-5.0).powf(2.0)
}

fn himmelblau(a: &Array1<f64>) -> f64 {
    if a.len() != 2 {
        panic!("input dimension must be 2.");
    }
    let x = a[0];
    let y = a[1];
    (x*x+y-11.0).powf(2.0) + (x+y*y-7.0).powf(2.0)
}

fn ackley(a: &Array1<f64>) -> f64 {
    if a.len() != 2 {
        panic!("input dimension must be 2.");
    }
    let x = a[0];
    let y = a[1];
    let pi = consts::PI;
    -20.0 * (-0.2 * (0.5 * (x*x + y*y)).sqrt()).exp() 
             - (0.5*(2.0*pi*x).cos()+(2.0*pi*y).cos()).exp() + 1.0_f64.exp() + 20.0
}

fn nelder_mead<T>(f: T, x_init: &Array1<f64>) -> (f64, Array1<f64>)
    where T: Fn(&Array1<f64>) -> f64
{
    // dimension
    let n = x_init.len();

    // constant
    let delta_e = 2.0;
    let delta_oc = 0.5;
    let delta_ic = -0.5;
    let gamma = 0.5;

    // initial simplex
    // https://github.com/scipy/scipy/blob/master/scipy/optimize/optimize.py
    let nonzdelt = 0.05;
    let zdelt = 0.00025;

    let mut y: Vec<(f64, Array1<f64>)> = Vec::new();
    y.push((f(&x_init), x_init.clone()));
    for k in 0..n {
        let mut y_k = x_init.clone();
        if y_k[k] != 0.0 {
            y_k[k] = (1.0 + nonzdelt) * y_k[k];
        } else {
            y_k[k] = zdelt;
        }
        y.push((f(&y_k), y_k));
    }

    for _ in 0..1000 {
        y.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());

        let f_best = y[0].0;

        // Reflect
        let mut _x_c = y[0].1.clone();
        for i in 1..n {
            _x_c += &y[i].1;
        }
        let x_c = &_x_c / (n as f64);

        let x_r = &x_c + &(&x_c - &y[n].1);
        let f_r = f(&x_r);
        if f_best <= f_r && f_r < y[n-1].0 {
            y[n] = (f_r, x_r);
            continue;
        }

        // Expand
        if f_r < f_best {
            let x_e = &x_c + &(delta_e * (&x_c - &y[n].1));
            let f_e = f(&x_e);
            if f_e < f_r {
                y[n] = (f_e, x_e);
                continue;
            } else {
                y[n] = (f_r, x_r);
                continue;
            }
        }

        // Outside Contraction
        if y[n-1].0 <= f_r && f_r < y[n].0 {
            let x_oc = &x_c + &(delta_oc * (&x_c - &y[n].1));
            let f_oc = f(&x_oc);
            if f_oc < f_r {
                y[n] = (f_oc, x_oc);
                continue;
            } else {
                y[n] = (f_r, x_r);
                continue;
            }
        }

        // Inside Contraction
        if y[n].0 <= f_r {
            let x_ic = &x_c + &(delta_ic * (&x_c - &y[n].1));
            let f_ic = f(&x_ic);
            if f_ic < y[n].0 {
                y[n] = (f_ic, x_ic);
                continue;
            }
        }

        // Shrink
        for i in 1..n+1 {
            let y_i = &y[0].1 + &(gamma * (&y[i].1 - &y[0].1));
            let f_i = f(&y_i);
            y[i] = (f_i, y_i);
        }
    }
    y.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
    (y[0].0, y[1].1.clone())
}

fn main() {
    let y = array![3.0, 3.0];
    let (min, v) = nelder_mead(&booth, &y);
    println!("a minimum of Booth: {} when (x, y) = ({})", min, v);
    let (min, v) = nelder_mead(&himmelblau, &y);
    println!("a minimum of Himmelblau: {} when (x, y) = ({})", min, v);
    let (min, v) = nelder_mead(&ackley, &y);
    println!("a minimum of Ackley: {} when (x, y) = ({})", min, v);
}

結果は、

a minimum of Booth: 0 when (x, y) = ([1.0000000000000004, 3])
a minimum of Himmelblau: 0 when (x, y) = ([3, 2])
a minimum of Ackley: 7.250114821582645 when (x, y) = ([2.987540988158525, 2.9937622971963957])

となり、Booth関数とHimmelblau関数については大域最小値(の一つ)を発見できました。
Ackley関数ではうまく動作しませんでしたが、これは関数の形状的に仕方がないでしょう。


以前のバージョン

この記事の以前のバージョンでは英語版Wikipediaに書いてあるアルゴリズムを参考にしたコードを載せていました。
当時ハマったところとその解決編を合わせて残しておきます。

extern crate ndarray;

use ndarray::*;

fn booth(x: &Array1<f64>) -> f64 {
    (x[0]+2.0*x[1]-7.0).powf(2.0) + (2.0*x[0]+x[1]-5.0).powf(2.0)
}

fn nelder_mead<T>(f: T) 
    where T: Fn(&Array1<f64>) -> f64
{
    // constant
    let alpha = 1.0;
    let gamma = 2.0;
    let rho = 0.5;
    let sigma = 0.5;

    // initial simplex
    let x0 = array![0.0, 0.0];
    let x1 = array![5.0, 0.0];
    let x2 = array![0.0, 5.0];

    let mut x = vec![(f(&x0), x0), (f(&x1), x1), (f(&x2), x2)];

    for _ in 0..100 {
        x.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());     // ハマりポイントその1
        println!("best value, vector: {}, {}", x[0].0, x[0].1);

        // centroid
        let xo = (&x[0].1 + &x[1].1) / 2.0;  // ハマりポイントその2a

        // Reflection
        let xr = &xo - &(alpha * (&xo - &x[2].1)); // ハマりポイントその2b
        let f_xr = f(&xr);
        if x[0].0 <= f_xr && f_xr < x[1].0 {
            x[2] = (f_xr, xr);
            continue;
        }

        // Expansion
        if f_xr < x[0].0 {
            let xe = &xo + &(gamma * (&xr - &xo));
            let f_xe = f(&xe);
            if f_xe < f_xr {
                x[2] = (f_xe, xe);
                continue;
            } else {
                x[2] = (f_xr, xr);
                continue;
            }
        }

        // Contraction
        let xc = &xo + &(rho * (&x[2].1 - &xo));
        let f_xc = f(&xc);
        if f_xc < x[2].0 {
            x[2] = (f_xc, xc);
            continue;
        }

        // Shrink
        for i in 1..3 {
            let xi = &x[0].1 + &(sigma * (&x[i].1 - &x[0].1));
            let f_xi = f(&xi);
            x[i] = (f_xi, xi);
        }
    }
}


fn main() {
    nelder_mead(&booth);
}

結果

cargo runすると

best value, vector: 9, [0, 5]
best value, vector: 9, [0, 5]
best value, vector: 5.9140625, [2.8125, 1.5625]
best value, vector: 1.30712890625, [1.328125, 2.265625]
(中略)
best value, vector: 0.1797569358295823, [1.3125408096108064, 2.7780022306583767]
best value, vector: 0.1797569358295823, [1.3125408096108064, 2.7780022306583767]
best value, vector: 0.1797569358295823, [1.3125408096108064, 2.7780022306583767]

となりました。Booth functionはf(1, 3)=0が最小値なので微妙です。

ハマりポイント(コード参照)

その1

Nelder-Mead法の実装には関数値でソートする必要があるのですが、f64型に対してはsort()を使うことができませんでした。
How to sort a Vec of floats? - help - The Rust Programming Language Forumを見てこう書きました。

追記: ordered-floatというのがあるらしいです。

その2

最初は2aを

let xo = (&x[0].1 + &x[1].1) / 2.0;

と書き、2bを

let xr = &xo - alpha * (&xo - &x[2].1);

と書いていたのですが、前者はコンパイルを通り、後者はコンパイルが通りませんでした。
エラーメッセージを読むと期待した通りの型になっていないっぽいのですが、どうしてなのかは分かりませんでした。

その後

let xr = xo.clone() - alpha * (xo.clone() - x[2].1.clone());

としたらコンパイルが通るようになりました。
「clone()すると複製コストがかかるような気がする(気のせいかもしれない)ので、cloneしない方法があるなら知りたいと思っています」と書いたところ、有識者からご指摘をいただき今のコードになりました。


Rustによるニュートン法の実装と可視化

Rustの勉強として、1変数関数の簡単なニュートン法を実装してgnuplotによる可視化を行いました。

可視化ライブラリ

GitHub - SiegeLord/RustGnuplot: A Rust library for drawing plots, powered by Gnuplot.

[Rust plot]で検索して一番上に出てきたこれを使いました。
正直gnuplotの表示は嫌いなのでmatplotlibやggplotのラッパーが欲しいです。

実装

cargo new newton --bin

で新しいプロジェクトを作り、Cargo.toml

[dependencies]
gnuplot = "0.0.24"

を追記しました。

 f(x):=x^2-2=0の解を求めるプログラムは以下の通りです。

extern crate gnuplot;
use gnuplot::*;

// 解きたい関数 f
fn f(x: f64) -> f64 {
    x*x - 2.0
}

// fの導関数
fn df(x: f64) -> f64 {
    2.0*x
}

fn main() {
    // 関数のプロットに使うx座標たち
    let mut t = vec![];
    for i in 0..100i32 {
        t.push(i as f64 * 0.1 - 5.0);
    }
    let mut fg = Figure::new();
    // axes2d()はmutable borrowを作るので後でshow()するには別スコープを作る必要がある
    {
        let axes = fg.axes2d();
        // x軸を表示
        axes.set_x_axis(true, &[]);
        // 表示範囲の指定
        axes.set_x_range(Fix(0.0), Fix(5.0));
        axes.set_y_range(Fix(-5.0), Fix(30.0));
        // fのプロット
        axes.lines(t.iter(), t.iter().map(|&x| f(x)), &[Color("red")]);

        // ニュートン法
        let mut x = 4.0; // 初期値
        while f(x).abs() > 1e-10 { // f(x)が十分小さくなるまで続ける
            // 関数値のプロット
            axes.points(&[x], &[f(x)], &[Color("blue"), PointSymbol('O')]); 
            // 関数値に縦線を引く
            axes.lines(&[x, x], &[-5.0, 30.0], &[Color("blue")]);
            // 接線のプロット
            axes.lines(t.iter(), t.iter().map(|&p| df(x)*(p-x)+f(x)), &[Color("black")]);
            // 値の更新
            x = x - f(x)/df(x);
        }
        println!("solution: {}", x);
    }
    fg.show();
}

cargo runの結果は

solution: 1.4142135623730951

で、以下のようなプロットが表示されました。
f:id:kivantium:20180210133822p:plain

Rustは初めて書くので変なところがあったら教えてください。

そのうち自動微分とか精度保証とかやっていきたいところ。

2018年1月第3,4週

2018年は隔週にします。

Webサイト

英語表現が実際に使われているか調べるサイト

科学論文だとSpringer Exemplar - Scientific Terms in Contextがおすすめ。

DNAシークエンスの値段グラフ

親の顔より見たグラフの元ネタはここっぽい。

RustでOSを書く

ツイート





論文

特定商取引法に定められた事項は請求により遅滞なく提供する