Search on the blog

2015年4月30日木曜日

scikit-learn(1) SVM

 SVMでirisを識別するサンプルコード。
各カーネルのパラメータはデフォルト値を使用。

clf.fit(X, y)で学習データ(特徴ベクトル, ラベル)を渡す。
clf.score(X, y)で識別率(※このサンプルでは学習データの識別率)を計算出来る。
from sklearn import svm, datasets

iris = datasets.load_iris()
X = iris.data
y = iris.target

# linear kernel
clf = svm.SVC(kernel='linear')
clf.fit(X, y)
print clf.score(X, y)

# RBF kernel
clf = svm.SVC(kernel='rbf')
clf.fit(X, y)
print clf.score(X, y)

# polynomial kernel
clf = svm.SVC(kernel='poly')
clf.fit(X, y)
print clf.score(X, y)

# sigmoid kernel
clf = svm.SVC(kernel='sigmoid')
clf.fit(X, y)
print clf.score(X, y)
実行結果は以下のとおり。

0.993333333333
0.986666666667
0.98
0.333333333333

2015年4月26日日曜日

おいら、レベル7に上がったっぺ


こんばんはごわす、おいら、レベル7に上がったっペ。

久々に問題を解きまくる休日を過ごした。
今日一日で七問解いた。やはり問題を解くのは楽しい。

今まであまり得意ではなかった”観察と分析を必要とするようなアドホックな問題”をいくつか解くことが出来て、数学のセンスが少し磨かれたような感覚を覚えた。

整数nが11で割れるかどうか

 整数nが11で割れるかどうかを考える。

10の冪乗の法11での剰余を見てみると

1 = 1 (mod 11)
10 = -1 (mod 11)
100 = 10 * 10 = 1 (mod 11)
1000 = 100 * 10 = -1 (mod 11)
10000 = 1000 * 10 = 1 (mod 11)
.....

のように1,-1,1,-1,...の繰り返しとなる。

よって、nの十進表記を"gfedcba"とするとnの法11での剰余は、

n = g - f + e - d + c - b + a (mod 11)

となる。

つまり、g - f + e - d + c - b + a = 0

であればnは11で割り切れる。

例)
3641 = - 3 + 6 - 4 + 1 = 0 (mod 11)より、3641は11で割り切れる.
13574 = 1 - 3 + 5 - 7 + 4 = 0 (mod 11)より、13574は11で割り切れる.

2015年4月21日火曜日

C++のconstの話

 少しトリッキーなconstの話。constキーワードを何処に付けるかによって
  1. constなデータを指すポインタ
  2. データを指すconstなポインタ
  3. constなデータを指すconstなポインタ
という違いが生じる。

百聞は一見に如かずということでコードを。
int main(int argc, char **argv) {

    const int x = 5;
    int y = 10;

    const int *p = &x;    // a pointer that points to a const int variable
    ++p;    // OK
    ++*p;   // NG

    int * const q = &y;  // a const pointer that points to an int variable
    ++q;    // NG
    ++*q;   // OK
    
    const int * const r = &x;   // a const pointer that points to a const int variable
    ++r;    // NG
    ++*r;   // NG

    return 0;
}

2015年4月19日日曜日

C++マルチスレッド入門

まえがき
先日のCode Jamで並列処理を行えばゴリ押しで解ける問題が出題された。
本番中ゴリ押し解を思いつくには思いついたのだが、C++でマルチスレッドの処理を書いたことが無くて、ごにょごにょやってるうちにタイムアップとなってしまった。

 せっかくなので、C++でマルチスレッドのプログラムの書き方を勉強してみた。後で見返した時にすぐ思い出すように簡単にまとめておく。
  1. マルチスレッド版Hello, World
  2. 子スレッドに引数を渡す
  3. mutexを使ったロック
  4. 子スレッドから結果を受け取る
  5. 子スレッドに後から引数を渡す

1. マルチスレッド版Hello, World
threadクラスのインスタンスを作成すると、スレッドが作成される。
スレッド化したい処理はコンストラクタで渡す。
コンストラクタに渡す処理は、関数でもファンクタでもラムダ式でもよい。
threadインスタンスのjoin()を呼ぶとスレッドの完了を待つ。
detachを呼ぶとスレッドの完了を待たずにmainは終了する。
#include <iostream>
#include <thread>

using namespace std;

int main(int argc, char **argv) {
    
    thread t1([]() { cout << "Hello, World!" << endl; });
    
    // do other tasks

    t1.join();       // wait t1 to finish
    //t1.detach();   // do not wait t1 to finish

    return 0;
}

2. 子スレッドに引数を渡す
スレッド化したい処理に引数がある場合は、threadのコンストラクタの第2引数以降で指定する。
引数を参照で渡したい場合は、refを使って参照ラッパーを渡す。
#include <iostream>
#include <thread>
#include <vector>

using namespace std;

void sayHello(vector<string> &names) {
    for (auto &x : names)
        cout << "Hello, "<< x << "." << endl;
}

int main(int argc, char **argv) {
    
    vector<string> names{ 
        "taro", 
        "hanako", 
        "jiro", 
        "mika"
    };

    thread t1(sayHello, ref(names));

    t1.join();

    return 0;
}
3. mutexを使ったロック
複数スレッドで共有するリソースを使う場合はmutexで排他制御を行う。
mutexクラスのlock、unlockメソッドを直接使うのは非推奨。RAIIでmutexの管理を行う標準ラッパークラスが提供されているのでそれを使う。
#include <iostream>
#include <thread>
#include <mutex>
#include <vector>
#include <chrono>

using namespace std;

mutex mtx_stdout;

void work(int t) {

    {
        lock_guard<mutex> guard(mtx_stdout);  // RAII
        cout << "work: " << t << " begin." << endl;
    }

    //....
    this_thread::sleep_for (chrono::seconds(3));
    
    {
        lock_guard<mutex> guard(mtx_stdout);  // RAII
        cout << "work: " << t << " end." << endl;
    }

}

int main(int argc, char **argv) {

    vector<thread> threads(10);

    for (int i = 0; i < 10; i++)
        threads[i] = thread(work, i+1);
    
    for (auto &x : threads)
        x.join();
    
    return 0;
}
より柔軟な制御を行いたい場合は、unique_lockを使うという選択肢もある。
void work(int t) {
    unique_lock<mutex> locker(mtx_stdout, defer_lock);

    locker.lock();
    cout << "work: " << t << " begin." << endl;
    locker.unlock();

    //....
    this_thread::sleep_for (chrono::seconds(3));
    
    locker.lock();
    cout << "work: " << t << " end." << endl;
    locker.unlock();
}

5. 子スレッドから結果を受け取る
子スレッドから処理結果を受け取りたい場合は、futureを使う。
asyncの第一引数にlaunch::asyncを指定することで別スレッドで処理を開始出来る。
第一引数にlaunch::deferredを指定した場合は、処理が結果取得時まで遅延される。(※別スレッドは作成されないことに注意)
#include <future>
#include <vector>
#include <iostream>

using namespace std;

bool isPrime(long long n) {
    if (n < 2)
        return false;
    
    for (long long i = 2; i * i <= n; i++)
        if (n % i == 0)
            return false;
    
    return true;           
}

int main(int argc, char **argv) {

    future<bool> fut = async(launch::async, isPrime, 1000000007);
    cout << fut.get() << endl;

    return 0;
}
6. 子スレッドに後から引数を渡す
スレッド実行時に引数の値が分かっておらず、後から非同期で子スレッドに引数を渡したい場合はpromiseを使う。
#include <future>
#include <vector>
#include <iostream>

using namespace std;

bool test(int x, future<int> &f) {
    return x < f.get();
}

int main(int argc, char **argv) {

    promise<int> p;
    future<int> f = p.get_future();

    future<bool> fut = async(launch::async, test, 10, ref(f));

    p.set_value(12);
    cout << fut.get() << endl;

    return 0;
}

スレッドに参照を渡したい場合はrefを使う

 C++でスレッドに参照を渡すときの話。

まずダメな例から。
#include <iostream>
#include <vector>
#include <future>

using namespace std;

int solve(vector<int> &x) {
    x[0] += 100;
    return x[0];
}

int main(int argc, char **argv) {
    
    vector<int> v{1,2,3};

    future<int> fut = async(launch::async, solve, v);

    cout << fut.get() << endl;
    cout << v[0] << endl;

    return 0;
}
上記のプログラムを実行すると、

101
1

と表示される。vector<int> v は参照渡しされているように見えるが、実際はされていない。

スレッドに参照渡しをする場合は、以下のようにrefを使って参照ラッパーを渡すようにしなければならない。
#include <iostream>
#include <vector>
#include <future>

using namespace std;

int solve(vector<int> &x) {
    x[0] += 100;
    return x[0];
}

int main(int argc, char **argv) {
    
    vector<int> v{1,2,3};

    future<int> fut = async(launch::async, solve, ref(v));

    cout << fut.get() << endl;
    cout << v[0] << endl;

    return 0;
}

実行結果は、
101
101
となる。

2015年4月18日土曜日

平方数を4で割った余り


平方数判定をする際にいくらか高速化できるテクニックを知ったのでメモしておく。

nが偶数の平方の場合は、
n = 2k * 2k = 4k^2
と表すことが出来る。よって4で割った余りは0。

nが奇数の平方の場合は、
n = (2k + 1) * (2k + 1) = 4k^2 + 4k + 1
と表すことが出来る。よって4で割った余りは1。

ということで、3とビットのandを取って1より大きくなった場合、その 数は平方数ではないことが分かる。

2015年4月10日金曜日

Lec 12 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分学12回目の講義をみた。このあたりは学生時代に研究室の先生方に丁寧に教えてもらったので、さすがに覚えている。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

講義の概要
  1. 勾配ベクトルの定義
  2. 勾配ベクトルは等高線と直交する
  3. 三次元曲面のある点における接平面の計算
    1. 勾配ベクトルを使った解法
    2. 一次の近似を使った解法
  4. 方向微分
  5. 勾配ベクトルは関数がもっとも増加する方向を向く

2015年4月5日日曜日

平方剰余の根

ある整数xに対して

x^2 = a (mod p)

となるようなaを「pの平方剰余」という。

以下ではaが与えられたときに

x^2 = a (mod p)

を満たすxを求めることを考える。ただし今回は簡単のためpが素数の場合のみを考える。

解の存在
ルジャンドル記号(a/p)を以下のように定義する。

a = 0 (mod p)ならば, (a/p) = 0
aが法pで平方剰余ならば、(a/p) = 1
aが法pで平方剰余でなければ、(a/p) = -1

pが奇素数の場合は、オイラーの基準により

(a/p) = a^{(p-1)/2}  (mod p)

が成り立つ。よって繰り返し二乗法を用いることで解の存在を簡単に判定出来る。
long long modpow(long long x, long long p, long long mod) {
    long long ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

// a is a quadratic residue modulo p?
int Legendre(long long a, long long p) {

    long long ret = modpow(a, (p-1)/2, p);

    if (ret == p-1)
        ret = -1;

    return ret;
}

解の公式
ルジャンドル記号が1の場合は、

x^2 = a (mod p)

を満たすxが存在する。

p=2の場合はxを全列挙すればよいため、以下ではpは奇素数とする。
pを8で割ったときの余りが3, 5, 7の場合は解の公式が存在する。

i) p = 3 (mod 8)、またはp = 7 (mod 8) の場合

x = ± a^{(p+1)/4} (mod p)

が解となる。

ii) p = 5 (mod 8)の場合

a^{(p-1)/4} = 1なら x = ± a^{(p+3)/8} (mod p)
a^{(p-1)/4} = -1ならx = ± 2^{(p-1)/4} a^{(p+3)/8} (mod p)

が解になる。
// calculate x such that x^2 = a (mod p)
// p is an odd prime that satisfies p = 3, 5, 7 (mod 8)
long long modsqrt(long long a, long long p) {
    
    // case i)
    if (p % 4 == 3)
        return modpow(a, (p+1)/4, p);

    // case ii)
    if (p % 8 == 5) {
        if (modpow(a, (p-1)/4, p) == 1)
            return modpow(a, (p+3)/8, p);
        else
            return modpow(2, (p-1)/4, p) * modpow(a, (p+3)/8, p) % p;
    }
        
    throw;
}
これで75%の場合は公式を使えば解けることが分かる。
しかし残りの25%の場合(p = 1 (mod 8)のとき)は公式が存在しないので、以下で示すアルゴリズムを使って計算するとよい。

非決定性アルゴリズム
Tonelli–Shanks algorithmを使うと、任意の奇素数pについて

x^2 = a (mod p)

を解くことが出来る。

以下に実装例を載せておく。
#include <iostream>
#include <random>

using namespace std;

long long modpow(long long x, long long p, long long mod) {
    long long ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

// a is a quadratic residue modulo p?
int Legendre(long long a, long long p) {

    long long ret = modpow(a, (p-1)/2, p);

    if (ret == p-1)
        ret = -1;

    return ret;
}

// Tonelli–Shanks algorithm
// calculate R such that R^2 = n (mod p)
// p is an odd prime
// n is a positive integer that satisfies (n/p) = 1
// See http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
long long modsqrt(long long n, long long p) {

    // step 1
    long long Q = p-1;
    long long S = 0;
    while (Q % 2 == 0) {
        ++S;
        Q /= 2;
    }

    if (S == 1) {
        return modpow(n, (p+1)/4, p);
    }

    // step 2
    default_random_engine generator;
    uniform_int_distribution<long long> distribution(0, p);
    
    long long z = 1;
    while (Legendre(z, p) != -1) {
        z = distribution(generator);
    }

    long long c = modpow(z, Q, p);
    
    // step 3
    long long R = modpow(n, (Q+1)/2, p);
    long long t = modpow(n, Q, p);
    long long M = S;

    // step 4
    while (t != 1) {
        long long i;
        long long t2 = t;
        for (i = 1; i < M; i++) {
            t2 = t2 * t2 % p;
            if (t2 == 1)
                break;
        }

        long long b = modpow(c, 1LL<<(M-i-1), p);
        R = R * b % p;
        t = (t * b % p) * b % p;
        c = b * b % p;
        M = i;
    }
    
    return R;
}

// test
bool isPrime(long long n) {
    if (n < 2)
        return false;
    
    for (long long i = 2; i * i <= n; i++)
        if (n % i == 0)
            return false;
    
    return true;           
}

#include <cassert>

int main(int argc, char **argv) {
    
    for (int i = 3; i < 10000; i++) {
        if (!isPrime(i))
            continue;

        for (int j = 0; j < i; j++) {
            if (Legendre(j, i) != 1)
                continue;

            long long x = modsqrt(j, i);
            assert (x * x % i == j);
        }

    }

    return 0;
}
参考URL
[1] Euler's criterion - Wikipedia, the free encyclopedia
[2] Computing square roots mod p
[3] Tonelli–Shanks algorithm - Wikipedia, the free encyclopedia

2015年4月3日金曜日

Lec 11 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分学の11回目の講義。
高校時代に暗記した公式の導出方法が出てくる。「ほう。あの公式はこうやって出すのか。」と目から鱗だった。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07
License: Creative Commons BY-NC-SA

講義概要
  1. 全微分と偏微分
    1. 微分の連鎖律
    2. 積の微分
    3. 商の微分
  2. 偏微分の座標変換(パラメータ変換)
2.ではデカルト座標(x, y)で表される関数の極座標(r, θ)パラメータに対する挙動を知りたい場合(もしくはその逆)どうするか?という例を紹介している。

2015年4月1日水曜日

リキッドレイアウトとレスポンシブデザイン

語句の意味
リキッドレイアウトとは、%を使ってウェブサイトをデザインすること。これにより、画面サイズが変わっても各エレメントが占める空間の割合が一定となる[1]。
 レスポンシブデザインとは、CSS Media Queries[2]を利用することで画面サイズに応じて異なるレイアウトを提供する手法のこと[1]。

デモ
まず普通のデザイン。各エレメントのサイズはpxで指定している。

sample1

画面サイズを小さくすると以下のように右側が隠れてしまう。


これをリキッドレイアウトに変更する。

sample2

すると、画面サイズを小さくしても右側が隠れない。各エレメントの幅の”割合”はディスプレイサイズによらず一定となる。


ちょっと詰まったところが、borderやpaddingの指定。divのwidthは%指定できてもborderはpx指定なのでサイズをピッタリ合わせるのが難しかった。
どうやらそういうときは、

box-sizing: border-box;

という指定をするらしい。この指定により、borderやpaddingが要素の幅に含まれるようになるらしい[3]。

リキッドレイアウトで画面サイズの問題は解決したように思えるが、さらにサイズを小さくすると以下のようになってしまう。


divエレメント内の文字が全部表示されていない。

これをレスポンシブデザインを使って解決してみる。

sample3

上のサンプルでは、CSS Media Queriesを使って画面サイズが500px以下になるレイアウトを変更するようにしている。これにより画面サイズを大幅に小さくした場合の対応が出来る。

今回の例では画面幅が500px以下になると、
  • 文字サイズ
  • メニューの段組
  • ヘッダーの文字色
を変更するようにしてみた。




参考ページ
[1] terminology - What is the exact difference between fluid and responsive design? - User Experience Stack Exchange
[2] CSS3 @media Rule
[3] box-sizing-CSS3リファレンス