ami_GS's diary

情報系大学院生の備忘録。ネットワークの勉強にハマっています。

matplotlibでboxplotの調べたことまとめ

はじめに

久しぶりです。研究やら授業でてんてこ舞いでした。

研究のためにboxplotで試行錯誤した結果のまとめを書こうと思います。



グループ化、及びbox内の塗りつぶし

普通にプロットすると、全部同じ色になってしまいます。
1つのグループに対し、複数のデータ群がある場合、同じ色だけで表現するのは非常に見難くなります。
なので、グループ毎にbox内を塗りつぶしてプロットしてみましょう。

今回は、3つの国にある、バナナ(黄色)ときゅうり(緑)の長さを例にしましょうかね、値は適当です。

from matplotlib import pyplot as plt
from scipy import stats

data = []
with open('./lengthData.txt', 'r') as f:
    for line in f.readlines():
        data.append(map(int, line.split())) # テキストデータから数値に変換してdataに格納

groupData = []
for i in range(3):
    tmpData = []
    for j in range(2):
        tmpData.append(data[i*2+j])  #(例えば) 3行2列の配列を作ります
    groupData.append(tmpData) #3サンプルそれぞれが2つのデータ群を持つ配列の完成

for i in range(3):
    #添字iでグループ指定、patch_artistはbox内の塗りつぶしです。
    #1グループの描画
    bp = plt.boxplot(groupData[i], positions=[2*i+1,2*i+2], widths=0.8, patch_artist=True) 
    bp['boxes'][0].set_facecolor('yellow') #黄色(バナナ)
    plt.setp(bp['medians'][0], color='red', linewidth=2) #メディアン
    plt.setp(bp['caps'][0], linewidth=1.5) #はこ髭の上限の線の太さ
    plt.setp(bp['caps'][1], linewidth=1.5) #はこ髭の下限の線の太さ
    plt.setp(bp['whiskers'][0], linewidth=2) #boxから伸びる上の点線の太さ
    plt.setp(bp['whiskers'][1], linewidth=2) #boxから伸びる下の点線の太さ

    bp['boxes'][1].set_facecolor('green') #緑(きゅうり)
    plt.setp(bp['medians'][1], color='red', linewidth=2)
    plt.setp(bp['caps'][2], linewidth=1.5) #同上(添字に注意!)
    plt.setp(bp['caps'][3], linewidth=1.5)
    plt.setp(bp['whiskers'][2], linewidth=2)
    plt.setp(bp['whiskers'][3], linewidth=2)

レジェンドを付ける

boxplotではデフォルトでレジェンドを付けることは出来ません。
(デフォルトのプロットの色が全部同じだから)
せっかく前節で色を付けたので、レジェンドを付けてみましょう。

banana, = plt.plot([1,1], 'y-', linewidth=10) #実際に線をプロットします。
cucumber, = plt.plot([1,1], 'g-', linewidth=10)

plt.legend((banana, cucumber), ('banana', 'cucumber')) #レジェンドの表示
banana.set_visible(False) #プロットした線を見えなくします。
cucumber.set_visible(False)

一度プロットしてからレジェンドを表示し、プロットした線を見えなくする、という手順です。




細かい設定

その他、プロットを見やすくするための設定です。

plt.ylim(0,40)
plt.yticks([i*20 for i in range(5)], ['0','10','20','30','40'], fontsize=17)
plt.ylabel('Length', fontsize=18)
plt.xticks([i for i in range(8)], ['','JP', 'JP', 'US', 'US', 'UK', 'UK', ''], fontsize=20)
plt.xlabel('Country', fontsize=18)
plt.grid()


以上の設定でプロットしたものが以下になります。
f:id:ami_GS:20140704103028j:plain



Pythonでしゃくとり法

はじめに

こんにちは、なかなか研究が軌道に乗らなくてツラいマンです。

今回はPythonでしゃくとり法を使う場面があったので載せていこうかなと。

しゃくとり法とは

与えられた配列の中から、ある条件を満たす最大の範囲を見つけるアルゴリズムです。

「左端(left)と右端(right)に配列の添字を入れ、rightをインクリメントしつつ、条件を満たさない要素を含んでしまった時に条件を満たさない要素が外にでるまでleftを進める」

という感じ?わかりにくいっすね

使用例

B: 細長いお菓子 - AtCoder Regular Contest 022 | AtCoder
を使います。

問題の解き方として、leftとrightの間に、同じ数字が入らないよう調整、その長さを取る。
という感じです

以下コード

N = int(raw_input())
A = map(int, raw_input().split())

left = 0 #左端初期値
ans = 1 #解初期値(0だとAの要素が1つだけの時にうまくいかない)
#rightはleftの1つ右からスタート
for right in range(1,N):
    #rightをインクリメントした後、同じ要素が含まれていた場合leftをインクリメント
    while A[right] in A[left:right]:
        left += 1
    ans = max(ans, right-left+1) #大きい方に更新
print ans


こんな感じです。
rightを進めてleftを縮める動きがしゃくとり虫っぽいからしゃくとり法なんですね、多分。

ちなみに

この解法だと最後のテストケースでTime Limite Exceededがでてしまい99点止まりです。
これを抜ける方法募集してます!(切実)


最後に

絶賛競技プログラミング練習中です。
アルゴリズムに苦手意識があるので覚えられる解法は全部覚える勢いで行っちゃいましょうね〜〜。

matplotlibの複数subplotをアニメーションで動かすヤツ

はじめに

こんにちは、GWですが周りの人の予定が合わなくて引きこもってるマンです。

題名の通り、matplotlibでの複数subplotをアニメーションで動かしていこうと思います。
研究で複数チャンネルの信号をリアルタイムで見たいと思い、書いてみたものです。

matplotlibとは?

簡単に言うとMATLABのplot系と同じような処理ができるようになるライブラリです。
データはnumpyの型をそのまま使えます。
少々速度に問題がありそうですが、MATLABが無かったり、サクッとプロットを出したい人には便利かも。

コード

matplotlibでアニメーションをする際には、matplotlib.animation.FuncAnimation()関数を使います。
引数にジェネレータ関数を取るので若干複雑(に見える)かも知れません。


コードは以下のとおり。

from matplotlib import animation
from matplotlib import pyplot as plt
import numpy as np

channelNum = 3       #チャンネル数
samplingrate = 100   #サンプリングレート
sample = 10          #1回のループで取るサンプル数
XMIN, XMAX = 0, 2    #x軸(時間を表す)
YMIN, YMAX = -1, 1   #y軸
INTERVAL = 1000/samplingrate*sample  #アニメーションさせる速さ(ms)
TIME = np.linspace(XMIN, XMAX, samplingrate*(XMAX-XMIN)) #x軸にサンプル点を置いておく

fig = plt.figure()
ax = []
plots = []
ydata = []

for i in range(channelNum):
    tmpax = fig.add_subplot(channelNum,1,i) #チャンネル数だけサブプロット作成
    tmpax.set_xlim(XMIN,XMAX)
    tmpax.set_ylim(YMIN, YMAX) #x,y軸設定
    ax.append(tmpax)
    ydata.append(np.zeros(0))
    plots.append(tmpax.plot(np.zeros(0), ydata[i])[0]) #チャンネルごとのプロット領域にデータをセット

def gen():
    while True:
        data = np.random.rand(sample, channelNum)*2-1  #チャンネルごとにsample分だけランダムな値を取る
        if ydata[0].shape[0] + len(data[:, 0]) >= len(TIME):
            initData() #プロットが右端まで達したらデータ初期化
        yield data

def initData():
    for i in range(channelNum):
        ydata[i] = np.zeros(0) #空にする

def updataData(data):
    for i in range(channelNum):
        ydata[i] = np.append(ydata[i], data[:, i])
        plots[i].set_data(TIME[:ydata[i].shape[0]], ydata[i])

    return plots  #データが格納された配列を返す

if __name__ == "__main__":
    #updateDataに第3引数のジェネレータ(gen)から返されたデータが引数として渡される
    #Macの場合はblit=Falseでないと動かない?
    ani = animation.FuncAnimation(fig, updataData, gen, blit=False, interval=INTERVAL)
    plt.show()

このような感じになります。

channelNum=6でプロットした物が以下になります。
f:id:ami_GS:20140504153901p:plain

アニメーションするのは1つで結構という方はchannelNumを1にしてください。
channelNumが3くらいであれば問題なく動くのですが、10くらいの大きさになると結構遅くなるので注意。

最後に

もしチャンネル数を多くした時の速度が気になるのであれば、tkinter, pyside等のGUIを使ってキャンバスに書き込んだほうが速いと思います。

今回はmatplotlibの紹介ということで許してちょんまげ!

Raspberry Piのアクセスポイント化 & ルータ化がうまくいかない時はこの通りにやればおk

はじめに

こんにちは。RPiのアクセスポイント化、ルータ化で詰まっていた情弱です。
今回、自分がはまっていた所を解決してくれた記事(英語)を見つけたので、それに沿って手順を書いて行きますね。

必要なステップ
  1. ワイヤレスアダプタ固定IPアドレス
  2. DHCPサーバのインストール及び構築
  3. アクセスポイントデーモン(hostapd)のインストール及び構築
  4. Wifiイーサネット間のルーティング設定


この例では、192.168.42.n のIPアドレスWifiネットワーク、192.168.1.nをイーサネット側のIPアドレスとして使います。

ワイヤレスアダプタ固定IPアドレス

"/etc/network/interfaces"を編集します。

sudo emacs /etc/network/interfaces

次に、
"allow-hotplug wlan0" 以外のwlan0、wpaに関する行をコメントアウトし、以下のように固定IPアドレスの設定をします。(/etc/network/interfacesの後半部分)

allow-hotplug wlan0
iface wlan0 inet static
address 192.168.42.1
netmask 255.255.255.0
gateway 192.168.1.1

#iface wlan0 inet manual
#wpa-roam /etc/wpa_supplicant/wpa_supplicant.conf
#iface default inet dhcp

DHCPサーバのインストール及び構築

DHCPサーバをインストールします

sudo apt-get install isc-dhcp-server

設定

sudo emacs /etc/dhcp/dhcpd.conf

"option domain-name"の部分(2行)をコメントアウトし、"authoritative"の部分をアンコメントします。

#option domain-name "example.org";
#option domain-name-servers ns1.example.org, ns2.example.org;

# If this DHCP server is the official DHCP server for the local
# network, the authoritative directive should be uncommented.
authoritative;

そして同ファイルの最後に以下のモノを加えてください

subnet 192.168.42.0 netmask 255.255.255.0 {
    range 192.168.42.10 192.168.42.50;
    option broadcast-address 192.168.42.255;
    option routers 192.168.42.1;
    default-lease-time 600;
    max-lease-time 7200;
    option domain-name "local";
    option domain-name-servers 8.8.8.8, 8.8.4.4;
    }

次に以下を設定します。

sudo emacs /etc/default/isc-dhcp-server

そして、INTERFACES=""  →  INTERFACES="wlan0"   と書き加えてください。

これを終えたら、次のコマンドでdhcpサーバを再起動します。

sudo service isc-dhcp-server restart

wifiドングルを指してからでないとdhcpサーバが起動しないので注意!

アクセスポイントデーモン(hostapd)のインストール及び構築

インストール

sudo apt-get install hostapd


以下の様にWifiネットワークの設定を書きます。

sudo emacs /etc/hostapd/hostapd.conf
interface=wlan0
driver=nl80211
#driver=rtl871xdrv  #最後に説明します
ssid=MyPi
hw_mode=g
channel=6
macaddr_acl=0
auth_algs=1
ignore_broadcast_ssid=0
wpa=2
wpa_passphrase=raspberry
wpa_key_mgmt=WPA-PSK
wpa_pairwise=TKIP
rsn_pairwise=CCMP

そして、上記の設定を/etc/default/hostapdから読み込むようにします。

sudo emacs /etc/default/hostapd
DAEMON_CONF="/etc/hostapd/hostapd.conf"

Wifiイーサネット間のルーティング設定

以下をアンコメントしてフォワーディング設定をします。

sudo emacs /etc/sysctl.conf
# Uncomment the next line to enable packet forwarding for IPv4
net.ipv4.ip_forward=1
echo 1 > /proc/sys/net/ipv4/ip_forward

以下のコマンドでRPiをルータにする設定をします。

sudo iptables -t nat -A POSTROUTING -o eth0 -j MASQUERADE
sudo iptables -A FORWARD -i eth0 -o wlan0 -m state --state RELATED,ESTABLISHED -j ACCEPT
sudo iptables -A FORWARD -i wlan0 -o eth0 -j ACCEPT
sudo iptables-save > /etc/iptables.ipv4.nat #ルーティングテーブルのセーブ

なお、hostapd.confを書き換えた場合はこれをもう一度行ってください


最後に、起動時にルーティングテーブルの読み込みを行うようにします。

sudo emacs /etc/network/interfaces

以下を最後の行に加えてください

pre-up iptables-restore < /etc/iptables.ipv4.nat

これでうまくいかない場合.....

“hostapd.conf”に"driver="が2行ありましたね、そこが肝になってきます。
それぞれが使っているUSBのWifiドングルにより、どちらを使うかが決まってきます。


どちらが使えるかは、以下のコマンドでわかります

sudo apt-get install iw
iw list

"nl80211 not found"
という出力が出た場合、"driver=nl80211"は使えません。"driver=rtl871xdrv"を使いましょう。


さらに、rtl871xdrvを使う場合は、使うべきhostapdも変わってきます。
以下のように設定するだけでOKです!

wget http://www.adafruit.com/downloads/adafruit_hostapd.zip 
unzip adafruit_hostapd.zip 
sudo mv /usr/sbin/hostapd /usr/sbin/hostapd.ORIG 
sudo mv hostapd /usr/sbin
sudo chmod 755 /usr/sbin/hostapd

これでアクセスポイント(SSID:MyPi)を見つけられるはずです。

以上。

日本語の記事では"driver=nl80211"を扱った記事がたくさんあったのですが、"driver=rtl871xdrv"を扱ったものが無くて全然構築が出来ずに困っていたわけです。

これで皆もRPiをアクセスポイントに出来る!やったね!!


Project euler #14 メモ化

3日に1記事挫折しそうやばいやばい

今回はProject eulerを題材にメモ化について書きます。

The following iterative sequence is defined for the set of positive integers:

n → n/2 (n is even)
n → 3n + 1 (n is odd)

Using the rule above and starting with 13, we generate the following sequence:

13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1
It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. 
Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.

Which starting number, under one million, produces the longest chain?

NOTE: Once the chain starts the terms are allowed to go above one million.

nが偶数:n = 3n + 1
nが奇数:n = n/2
を繰り返して、nが1になるまでのステップが最大になるnは何か、(n <= 1000000)
という問題。

単純に書くと

max = 0
n = 0
upper = 1000000

def length(n):
    if n == 1:
        return 1
    elif n%2 == 0:
        return 1 + length(n >> 1) #n/2のシフト演算版
    else:
        return 1 + length(3*n+1)

for i in range(1, upper):
    step = length(i)
    if step > max:
        max = step
        n = i

こんなバカ正直に書くと、僕の環境では計算終了まで30秒弱もかかります。
おそスギィ!!!!




メモ化手法

問題文に

13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1

とあるように、
n=13なら10ステップ、n=40なら9ステップ、n=20なら・・・と解が続いていきますね、


例としてあげたコードのように、1から総当たりで解を探していく場合、
n=5ならば6ステップという解が出ているのに、n=10の場合に、

10 → 5 → 16 → 8 → 4 → 2 → 1

のように計算し続けてステップ数を数えるのは非効率です。
そこで、n=5ならば6ステップという解をどこかに補完しておけば、

10 → 5

のように、1回の計算だけで「n=10ならば7ステップ」という解が出せます。




実装

(この場合は)実装は単純です。

memo = {}
max = 0
n = 0
upper = 1000000

def length(n):
    if n in memo:
        return memo[n]#メモにnの値が記録されていたらそのvalueを返す
    if n == 1:
        return 1
    elif n%2 == 0:
        return 1 + length(n >> 1)
    else:
        return 1 + length(3*n+1)

for i in range(1, upper):
    step = length(i)
    if step > max:
        max = step
        n = i
    memo[i] = step#メモにステップ数を記録

これだけ!

これを実行すると3秒切ります。

これだけで1/10も速くなるなんてすごい!

ー終ー

Raspberry Pi初期設定(ブート、ssh設定、固定IP)

卒論提出完了しました。これで自分の勉強に集中できるぞ!(3週間後に発表あり)


今回は、Raspberry Piの初期設定について書こうかと思います。
研究室の先生が僕に貸してくれてありがたやありがたや・・・


必要な物

作業環境Mac OS X
Raspberry Pi (Type B)
SD カード (今回はclass4 8GBの物を仕様)
電源用のアダプタ&マイクロUSBケーブル
HDMIケーブル&ディスプレイ
LANケーブル
キーボード
GUI操作をするのであればマウス


※microSDカードをSDカードの大きさにできるアダプタ?を使うと起動時に認識されないので注意



全部装着した物がこちら
f:id:ami_GS:20140125010143j:plain
電源が入ると左下のLEDが光り輝きます。

続きを読む

Pythonにおける*及び**の役割

3日に1回書く宣言してたけど人生稀に見る多忙さ故に挫折しそう。

なので簡単な事かきます。



色々な人のソース見てると、たま〜に、

def func(*foo, **bar):

このような物を見かけると思います。
Cでつまずいた人だと、ポインタだと思ってびっくりしちゃうかもしれません(びっくりしました)

結論から言うとポインタではありません。安心してください。

Pythonにおけるアスタリスクは可変長引数の役割をします。
引数の数に制限を欠けたくないときに便利ですね。

まず、*の場合(一つ)

def func(foo, bar, *args):
    print foo+bar

    for arg in args:
        print arg ,

func(1,2,"a","b","c")

とすると、出力は

3
a b c

となります。

foo及びbarに入る引数以外はタプルになります!便利ですね。
ちなみに入力が0個になっても動きます。


次に、**の場合(二つ)

def func(foo, bar, **kwargs):
    print foo+bar
    print kwargs

func(1,2,a="10",b="20",c="30")

とすると、出力は

3
{'a': '10', 'c': '20', 'b': '30'}

となります。辞書型ですね。


はい、見たところ確実にポインタではないようです。