某での某某某

僕は怒っています。
就活のコーディングテストでよく使われる某某サイトの入力がゴミウンチ以下だからです。
普段競プロでしかほとんどプログラミングしなかったり、あんまり文字列扱わなかったりすると変な入力が来たらウンチを漏らしてしまうものです。

そこで某の注意点を少しだけまとめます。
主にc++想定です。Pythonも触れるかも。(あんまり触れませんでした)

言語が途中で変えられない

見ての通りです。しかしオンラインジャッジということもありTLE(時間オーバー)はあります。
問題をみて言語変えられないのは僕的に辛い。
某社のテストの時は、pythonを選んでしまいめちゃくちゃ最適化したものの、もう原理的に時間内に終わらないやろってレベルで無理でした。

結論
C、C++は正義、速度はパワー
でも難しい処理を普段Pythonに任せていると辛いものがある。

入力がcharだしクソ(僕が受けた範囲では)

これ僕のような人間からしたら面食らいます。(競プロの標準入力は優しい)
数字をどうこうしようってアルゴリズムを組めと言われたので、入力が数字だと思い込んで実装した後に気付く・・・ということもありました。
Pythonであればテキトーにやるだけで変換できるので特に問題はありません。
しかし、改行のエスケープ文字とかが文末に入っててややこしいことになったりします。(入力に含まれるのでそのままforとかで回すと色々やばい)
正直こんなことで頭を使いたくないのでcharで受け取ってintにキャストするやつを置いておきます。(テスト中に焦りながら調べた)

#include <stdlib.h> /* strtol */
#include <limits.h> /* INT_MAX, INT_MAX */

int main(){
char *endp;
char *str;
long r;

r = strtol(input, &endptr, 10);

if (input == endptr) {
   //文字列の処理;
}

if (r > INT_MAX || r < INT_MIN) {
   //オーバーフローの処理;
}

int val = (int)r; //intにキャスト
}

inputが標準入力で最終的にはvalが得たい値です。
atoiは？？と思うかもしれませんが、あれは文字列をかますと0を返すはずなので入力が0の場合と分解できません。
後エラー時の挙動が環境で違うとかで危ないらしい。

少なくとも上のソースで僕は通したことがあるので大丈夫なはず。

感想を書かせてくる

これもわりとありがちなんですよね。
書かすのはいいんですよ、書かすのは。
ただ実装でカツカツのレベルの問題にぶち当たった時がクソ

「やべっ、時間足りない〜〜〜」

って時にあっ、感想(実装の注意点とか苦労した点とか書かせてくる)書かなきゃとなるわけです。
まぁこれに関しては初見でもどうにかできるとは思いますが、最初から計画を立てて置いたほうが良さそうですね。
事前説明にも多分書いてあるけど

終わりに

とまぁここまで書きましたが僕はもうこれといって頑張る必要もないのでみなさん頑張ってください。

データ解析諸々まとめ(R編)

例にもよって自分用、Rのコードでまとめていきます。

多分同じ内容でPython編も書きます。

以下目次

・はじめに

この記事ではパッケージの導入から基本操作をそれぞれまとめていく。
まだまだ触ったことがない手法が多いので使ったらどんどん追記していくスタイルにしておく。

あと、まだ使ってないけど使いたいやつは空でも章だけ作ってます。

~表記上の決まりごと~
train:訓練データ
test:テストデータ
y:目的変数

~環境(2018/7/11時点)~
・MacbookPrp HighSierra :  Rstudio ver.1.1.453 , R ver.3.5
・Windows10 : Rstudio ver.1.1.453 , R ver.3.5

あんまり解析手法の理論を理解していないのでヤバいこと書いちゃうかも>

・便利パッケージの導入と使い方

まずは補助的なパッケージの導入について

基本的に導入に詰まったところはRとRstudio のアップデートで解決した。

・R自体の更新方法
Rが3.5.0へメジャーアップデートしたので簡単アップデート - 清水誠メモ
・Rstudio の更新方法
HelpからCheck for Updates

・dplyr

フィルタリングとかグルーピングとかデータフレームの操作が早いやつ。
C++で書かれてるかららしい。

~導入手順~

1. 各種ソフトアップデート確認
2. R(Rstudio)のコンソールでinstall.packages("tidyverse")
3. library(tidyverse)で確認

~基本操作~
ここに書くと膨大になってしまうので参考にした記事のリンクをはる。
dplyrを使いこなす！基礎編

行ごとに操作するのは簡単だが、列をある条件で動かしたりするのは厳しいっぽい。(多分、誰か知ってたら教えて)
例えば定数のカラムを消したい場合は

train_var <- train%>%dplyr::summarise_all(funs(var))#分散をそれぞれのカラムで計算し、代入
dropcol<-c()#落とすカラム名を入れていくやつ
for(i in 1:length(train)){
  if(dtrain_var[1,i]==0){
    dropcol<-append(dropcol,names(train)[i])#分散が０ならカラム名を追加
  }
}
train<-train[,-dropcol]#落とすやつ以外を訓練データに回す。

あと、%>%の記法もこいつらの中に入ってるぽい？
これは右に左のやつを引数として渡すという意味

scaled_data<-data%>%scale()

上のようにdplyr::を使わなくても良い(例はあるデータを正規化するコマンド)
可読性は高いか低いかわからないけどユーザーは割と見る。

・caret

ハイパーパラメータのグリッドサーチしてくれたりする。
対応してないパッケージもある(LightGBMとか・・・)

~導入手順~

1. 各種ソフトアップデート確認
2. R(Rstudio)のコンソールでinstall.packages("caret")
3. library(caret)で確認

~基本操作~
これも膨大になるので参考記事を貼って主要部分のみ説明する。
Rによる機械学習：caretパッケージの使い方 | Logics of Blue

今の僕が詰まったのはtuneGrid
環境的な問題かもしれないけど、グリッドするパラメータを全指定する必要がある。

tuneGrid = expand.grid(mtry=1:20,splitrule="variance",min.node.size=5)

上はまた後ほど紹介するランダムフォレストのパッケージrangerのもの。(同様の内容を下でまた書きそう)
mtryは１から２０の間でグリッドサーチ、splitrule、min.node.sizeは固定(最適化しない)
中身に適当なこと書いてエラーを見ると何を指定すべきか見れる。

~例文~
一応rangerでの例を挙げておきます。

modelRanger <- train(
  formula = y ~., 
  data = train, 
  method = "ranger", 
  tuneLength = 4,
  preProcess = c('center', 'scale'),
  num.trees = 500,
  tuneGrid = expand.grid(mtry=1:20,splitrule="variance",min.node.size=5),
  trControl = trainControl(method = "cv"),
  importance = 'impurity'
)

パラメータ名	デフォルト値	説明	推奨
method	-	解析パッケージ指定	-
tuneLength	-	各パラメータを何通りずつ試すか。	-
preProcess	-	解析データの加工、この場合は正規化	-
tuneGrid	-	グリッドサーチするパラメータと範囲をあらわに指定	-
trControl	-	どのように解析するか、この場合はクロスバリデーション	-

・doParallel

caretを使う場合にはセットで読み込んでおいたほうが良さそうな並列化のパッケージ。
現段階ではどの範囲で並列してくれるのかわかってません。(caret以外で)

~導入手順~

1. 各種ソフトアップデート確認
2. R(Rstudio)のコンソールでinstall.packages("doParallel")
3. library(caret)で確認

~基本操作~
caretを使う前に以下を実行する

cl <- makePSOCKcluster(detectCores())#何コアで並列するか
registerDoParallel(cl)#並列するコアの登録？

detectCores()で使っているパソコンの最大コア数を取得できる。
上の例は実行環境での最大並列

また、原因はよくわかっていないが同じコードでもMac(4コア)では動かずWindows(8コア)では動いたことがあった。
初期化がどうこうとかよくわからないエラーが出るが、スペックの問題と思われる。(検索してもあんまりわからなかった)

・解析用パッケージの導入と基本操作

解析を実行するパッケージの導入方法と基本操作をまとめる。
この場合も各種アップデートは必須。あとは64bitじゃないパソコンはドブに捨てましょう。

・線形回帰

実はあんまりやったことないので使ったら書きます。

・主成分解析(PCA)

よくPCA(Principle Component Analysis)って略されるやつですね。
分散が大きい特徴量であるほどデータの本質を捉えているはず！という仮定から生まれた手法です。
直感的にもまぁそんな気がするなぁという。
そこで特徴量の分散が最大になる軸でデータを再構成し直す解析手法です。この辺の説明はネットに無限に転がっているので検索しよう。

一番参考にしたサイト
Practical Guide to Principal Component Analysis (PCA) in R & Python

~導入手順~
デフォルトで入ってます。

~基本操作~

result.pca <- prcomp(train,scale=T)

ほとんど上のコマンドで尽きています。
trainは解析したいデータで、scaleは正規化するかしないか。

~詰まったところ~
定数列があった場合スケールができないのでエラーが出る。
→上に書いた処理(分散が0の列を落とす)などを挟んでどうにかする。

~実行後~

#old.op <- options(max.print=999999)#特徴量がバカ多い場合は表示が溢れるのでこれをやる
sink("component.txt")#標準出力先を切り替える
round(result.pca$rotation, 3)#主成分を出力
sink("summary.txt")
summary(result.pca)#主成分が何パーセントデータを説明しているかがわかる。
sink()

これのCumulative Proportionがそう

~その他応用~
主成分解析→回帰に使う場合は目的変数を外して主成分解析、９割以上(一般的に何パーがいいのかは知らない)の主成分を拾ってきて回帰させる。

train <- data.frame(result.pca$x)
train<-train[,1:m]#mまでの主成分を拾ってくる(結果を見て変える)

これで適当な回帰にぶっこむ

・決定木系

決定木系の解析がkaggleとかだと一番使われてそうと思っている。
(決定木をスッキリ説明する語彙力がない)
理論の気持ちとかも充実させていきたいから勉強したらその辺も追記したい。

以下にわかりやすい図とかが乗っている気がします。
KaggleのHousePrices問題を決定木系のアンサンブルで解く - mirandora.commirandora.com

あと、決定木系は線形分離問題に弱いらしい。
そのためかXGBには線形用っぽい関数があります。

・決定木

実は使ったことがない、使ったら書きます。

・ランダムフォレスト(randomForest、ranger)

アンサンブル学習の一つ。弱学習機を〜〜って説明されます。
簡単に説明すると、多数決みたいなもんです。アンサンブル学習は全部多数決(アンサンブル平均をとる)ですけどね。
この場合は多数決取る人たち(決定木)は変わりません。変わる奴がXGBとかLGBMとかですね。

~導入手順(randomForest)~

1. 各種ソフトアップデート確認
2. R(Rstudio)のコンソールでinstall.packages("randomForest")
3. library(randomForest)で確認

基本的に上のやつでいいとは思いますが、今はもっと高速な奴があります。
rangerって奴です。データが小さい場合は上のものでいいですが、特徴量が多い場合やデータが膨大な場合はこれが早いです。
あとcaretにも対応しているので使わない手はない(？)

速度比較
[1508.04409] ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R

~導入手順(ranger)~

1. 各種ソフトアップデート確認
2. R(Rstudio)のコンソールでinstall.packages("ranger")
3. library(ranger)で確認

説明書(デフォルト値などはここで確認、しかしnum.treesのデフォルト値がよくわからない)
https://cran.r-project.org/web/packages/ranger/ranger.pdf

~基本操作~
・最小構成

ranger=ranger(y ~. , data=train.data)

分類も回帰もひとまずはこれで実行可能。分類の場合一応fcator(y)をかけといた方がいいかも？

・そこそこじっくり使うくらいの構成(回帰を想定した数値)

ranger = ranger(
    formula = y ~ .,
    data = train,
    num.trees = 100,
    mtry = 5,
    splitrule="variance",min.node.size=5,
    write.forest = TRUE,
    importance = 'impurity'
  )

これらのパラメータは重要そうなので説明します。(説明できないところはわかったら追記)

パラメータ名	デフォルト値	説明	推奨
num.trees	-	決定木の数	-
mtry	特徴量のルート個数	1つの決定木がいくつの変数を持つか	デフォルトの√n個らしい
splitrule	"variance"	-	-
min.node.size	分類は1回帰は5	-	-
write.forest	True	オブジェクトの保存をするかしないか	Trueがいいと思います。
importance	-	変数の重要度をどう表すか？	-

赤くしたパラメータがcaretでグリッドサーチできる(しなければならない)ものです。

・caretでのグリッドサーチ

modelRanger <- train(
  formula = y ~., 
  data = train, 
  method = "ranger", 
  tuneLength = 4,
  preProcess = c('center', 'scale'),
  num.trees = 500,
  tuneGrid = expand.grid(mtry=1:20,splitrule="variance",min.node.size=5),
  trControl = trainControl(method = "cv"),
  importance = 'impurity'
)

上で同じコードを例に出していますが、なるべく章ごとに閉じた作りにもしたいのでもう一回貼ります。

パラメータ名	デフォルト値	説明	推奨
method	-	解析パッケージ指定、今回はranger	-
tuneLength	-	各パラメータを何通りずつ試すか。	-
preProcess	-	解析データの加工、この場合は正規化	-
tuneGrid	-	グリッドサーチするパラメータと範囲をあらわに指定	-
trControl	-	どのように解析するか、この場合はクロスバリデーション	-

ランダムフォレストは結構使うので充実させたいね。

・Xgboost

まだ使ってないので、詳しくは使ったら。
導入と軽い調べ物だけしたので

~導入手順(Xgboost)~

1. 各種ソフトアップデート確認
2. R(Rstudio)のコンソールでinstall.packages("xgboost")
3. library(xgboost)で確認

説明書
https://cran.r-project.org/web/packages/xgboost/xgboost.pdf

参考(パラメータなどについて)
xgboost のパラメータ - puyokwの日記

・LightGBM

XGBに対抗してマイクロソフトが作ったとか言われてるやつ。
基本はXGBと一緒らしい。そのライト版かつ正確な(と主張している)解析手法。
問題設定にもよりそうだけどkaggleでXGB と共にめちゃくちゃ流行ってるみたい。

残念ながらcaretには対応してないっぽい。

~導入手順(LightGBM)~
結構めんどい。
LightGBM/R-package at master · Microsoft/LightGBM · GitHub
一応ここに全部書いてあるので確認してください。あとMacでやったのでMac想定です。
gitとかもう入ってるよ〜って人はそこは飛ばして次に行ってください。

1.各種アップデート(R,Rstudio,Homebrewなど)
2.gitコマンドをターミナルで使えるようにする。

brew update
brew install git

3.Cmakeなどを使えるようにする。
gccがインストールされていないと多分動きません。

brew install gcc

この時バージョンも確認しておきましょう。僕は脳筋なのでgcc-7,gcc-8と打っていってコマンドが認識されたバージョンとしました。

brew install cmake

思ったよりインストールが長かった。待とう
4.公式のコマンドを入力

git clone --recursive https://github.com/Microsoft/LightGBM
cd LightGBM/R-package
# export CXX=g++-7 CC=gcc-7  # for macOS (replace 7 with version of gcc installed on your machine)
R CMD INSTALL --build . --no-multiarch

これでエラーでなければ多分大丈夫
5.Rのコンソールでlibrary(lightgbm)を確認

各種パラメータ
LightGBM/Parameters.rst at master · Microsoft/LightGBM · GitHub

~基本操作~
正直日本語でかつRで使ってる人全然いない(記事を書いてくれていない)ので、細かいところまでわかってないです。

・最小構成(最小もよくわかってないからとりあえず空)

empty

追記するので、それまではこの記事みてください、下の方にRで簡単なirisの分類をしてます。
LightGBM ハンズオン - もう一つのGradient Boostingライブラリ

・僕が今使ってるやつ(こっちは回帰)

#訓練用と評価用に分割
set.seed(123)
inTrain <- createDataPartition(y, p=.9, list = F)#caretの関数、均等な分割をしてくれるらしい
tr <- train[inTrain,]
va <- train[-inTrain,]
tr_ta <- y[inTrain]
va_ta <- y[-inTrain]
#lgb用のデータセット作成
lgb.train = lgb.Dataset(data.matrix(tr), label = tr_ta)
lgb.valid = lgb.Dataset(data.matrix(va), label = va_ta)

#ここから解析
params.lgb = list(
  objective = "regression"
  , metric = "rmse"
  , min_data_in_leaf = 5
  , min_sum_hessian_in_leaf = 100
  , feature_fraction = .5
  , bagging_fraction = .5
  , bagging_freq = 4
)

lgb.model <- lgb.train(
  params = params.lgb
  , data = lgb.train
  , valids = list(val = lgb.valid)
  , learning_rate = 0.003
  , num_leaves = 180
  , num_threads = 2
  , max_depth = 8
  , reg_alpha = .3
  , reg_lambda = .1
  , min_child_weight = 33
  , zero_as_missing = T
  , nrounds = 5000
  , early_stopping_rounds = 600
  , eval_freq = 50
)

実行すると環境によるかもしれませんがWarningがアホほど出てきますが、無視してオーケーです。
ちなみにこんなやつ
[LightGBM] [Warning] No further splits with positive gain, best gain: -inf

ハイパーパラメータの推奨値や意味はある程度まとまったら追記します。

パラメータ名	デフォルト値	説明	推奨
-	-	-	-
-	-	-	-
-	-	-	-
-	-	-	-
-	-	-	-

・SVM

使ったら書きます。

getch・getcheに関するめも(Mac/Unix用) 追記4/12(キーボード入力まわりについて)

とあるプログラムをCで書いてたら

#include <conio.h>
getch()

の部分でconio.hがないよと怒られる。

調べるとmac(というかunixには実装されてない模様)

他の手立てもあるらしいが、あんまり難しいことを考えたくないため、どうにかできないか検索

nonberanai.blog.fc2.com

僕「なーんだ、もうあるじゃんw簡単簡単w」

と言っていたのも束の間、includeフォルダにgconio.hを移動できない

su使って移動させればええやろ！と思ったが、それでもPermission deniedされる。

（ちなみにsuを有効化するのも１手順必要だった
Mac でルートユーザを有効にする方法やルートパスワードを変更する方法 - Apple サポート
）

おかしいなと思いつつ調べるとusr以下は最近のアップデートでルートユーザーでもアクセスできなくなっている模様

rootlessとかいうやつで。

やばくないですか？(やばくはない)

というわけで解除しました

qiita.com

これの通りにやればOKです。

ただ、rootlessはマルウェアからの防衛的な意味もあるようなので

作業が終わった同じ手順でcsrutil disableの部分を

csrutil enable

というコマンドにして有効化しておきましょう。

終わり。

追記4/12
・kbhitについて

#include <stdio.h>
#include <conio.h>

int main(){
	int key;
	int flag = 1;

	while(flag){
		if(_kbhit()){
			処理１
			flag = 0;
		}
	}
}

キーボード入力がヒットする→キーボード入力があると処理をするような関数

これもMacにはデフォルトでないみたいです、上のヘッダファイルにもありませんでした。

そこでここ参照

tricky-code.net

解決

・矢印キーの入力について
特殊入力はなんか色々特殊らしい、安直に調べたやつではうまく認識できなかった

Console Operation for RT-Series

結局ここで解決した。以下例文(unix用)インデントおかしいけど許して

#include <stdio.h>
#include <gconio.h>
#include <termios.h>
#include <unistd.h>
#include <fcntl.h>


int kbhit(void)
{
    ~略~
    }

    return 0;
}

int main(){

	switch(getche())
   {
    case 0x1b:
    switch(getche())
      {
			case 0x5b:
			switch(getche()){
			        case 0x41: printf("↑\n"); break;
		                case 0x42: printf("↓\n"); break;
		                case 0x44: printf("←\n"); break;
		                case 0x43: printf("→\n"); break;
			}
      }
    break;
   }

return 0;
}

ボルツマンマシン

この前のこれ

の記事を書いていこうと思います。

はじめに

動画のこいつは何者なのか、とかいう質問が出てくると思いますが最初に答えます。

人類に限りなく近いゴリラです。

本題に入りましょう、今回のモデルは情報分野だと「ボルツマンマシン」とか「ホップフィールドネットワーク」とか言われてるやつです。

物理科の民は「え、ただのイジングモデルじゃん・・・。」ってなるやつです。

この辺の説明は後でやります。

流れとしては、説明書いてコードはって結果はって終わり！ってなると思います。

各モデルの説明

上で出てきたモデルの説明をしていきます。

物理の話から情報分野に下る流れでいきます。

・イジングモデル(Ising-model)

統計物理とかやった人はわかると思いますが、磁性体のモデルの中で一番(?)簡単なやつです。

まず、棒磁石を思い浮かべてください。おそらく大勢の人は赤と青に塗られ、それぞれにN、Sと書いてあるやつを思い浮かべると思います。

こいつらを半分に切ったらどうなるでしょう？

答えは「長さが半分になった棒磁石が二つできる」です。知ってる人も知らない人もいると思いますが・・・。

では、それぞれもう半分にしたらどうなるでしょう？

もう簡単だと思いますが「さらに半分になった棒磁石が四つできる」ということになります。

これを繰り返していくと、これ以上分割できない「磁石の最小単位」にたどり着くことがわかると思います。

これが「スピン」です。

このスピンは矢印に例えられ、矢印がみんなで同じ方向を向いていると全体としてN極だとかS極だとかが生まれるわけです。

イジングモデルというやつもこれをモデル化したもので、矢印は上か下かのどちらかしか向きません。

全部上に向いていればその方向がN極で、逆はS極ということになります。

大体こんなイメージです。全部上向いてますね

この時矢印(スピン)がどこを向きたいかと言うものを指定しなければいけません。

それがハミルトニアン(エネルギー関数)です。今回の場合

こう書けます。本来は量子力学の波動関数の重なりがどうとか議論して積分したりして得られるものですが、直感的にも良さそうなことがすぐにわかります。

またSはベクトルで書いていますがイジングモデルにおいては±１しかとりません(上か下か)、あとsummationはあらゆるijの組で取ることにしています。

例えばJijが全部＋１の場合を考えましょう。しかし、どういう状態が「答え」でしょうか？

現実世界というのは基本的に万物が全体としてエネルギー安定な方向へ向かっています。高いところから物を落とすと下に落ちるのもそうです。

何が言いたいかというと、このエネルギー関数は現実では常に最小化される状態(向き)になるということです。(厳密には自由エネルギーが最小化)

全てがJij=＋１となる時は簡単で全てが上、または下に揃えばエネルギーが最小になることがわかると思います。

これはスピンが何個ある時でも同じです。(N＞２)　磁石は全体として同じ方向を向いているので、こういったJijによって振る舞いが決まっていそう、ということも納得できると思います。

例ではJijが全て等しい時を考えましたが、Jijが色々な値を取るとそれぞれの矢印がもはやどういった向きをとれば良いのかわからないと思います。実際この場合は解析的に解くことはできません。

・ホップフィールドネットワーク(Hopfield-network)/ボルツマンマシン(Boltzmann-machine)

まずはホップフィールド ネットワークについてです。これはもうイジングモデルそのものです「スピンをJijと言う結合定数で結んでいる」と考えるとなんとなくネットワーク感も出てくるのではないでしょうか。具体的には下図のようになっています。

Jijの値が有限＝繋がっている(繋がってないところは０)と考えるとイジングモデルはこのようになり、情報分野ではホップフィールド ネットワークと呼ばれます。

また、CNNなどは順方向に情報が伝わっていき、ニューロンが次々と発火していくような向きを持った構造になっています。しかしこのモデルは向きを持たない無向グラフの一種です。

次にボルツマンマシンについてですが、ホップフィールド とほとんど等価です。何が違うのかと言うと、「温度」の概念が取り込まれていることです。

温度を取り込む利点として、記憶を想起させる時に便利とかはあると思いますが、他にはあるんでしょうかねえ〜〜勉強してないのでよくわかりません。

・結局何が違うのか

結局一緒なんじゃん、となると思います。はっきり言うと

結局一緒

です。本質的には何も変わりません。

じゃあ情報と物理で何が違うかと言うと「問題設定」です。

物理では、物質構造や現象から「Jijを予め決めてスピンの向きを決める」のに対して

情報では基本的に「スピンの向きからJijの組みを推定する」と言う問題を解きます。

逆問題の関係になっているわけですね。

例えば今回の話だと「２値画像(１ピクセル＝１スピン)からJijを推定(学習)」し、それを今度は物理の問題として「与えられたJijの組みに対してのスピンの向きを決める(学習画像を想起する)」と言う二段階になっています。

学習・アルゴリズム

ここまででイメージがついたと思うので(思うので)実装に入ります。

正直実装にあたって触れていない重要なことがあるんですが、１つの記事に書ききれないのでいくつか天下り的に紹介します。気になったら調べてください。

学習方法

ホップフィールド における二値画像の学習方法は古くから知られていて「Hebb則」ってやつを使います。

Jij:i番目とj番目のピクセルの結合定数
N:総ピクセル数
p:覚えさせたいp番目の画像
S:2値をとる学習させたい画像のピクセル

こんな感じです。今回の設定だと与えた画像のセットから決定論的に決まるので推定とかではないですが、立派な学習です。

また、Jiiなど添え字が同じやつは０にします。(自分自身への結合は０ということ)

想起方法

「メトロポリス・ヘイスティングス法」を使います。

これは簡単に言うと「乱数で１つのピクセルの向きを決めて、それがエネルギーを下げる方向ならば採用する」というのを繰り返すモンテカルロ法の１種です。ボルツマンマシンの文脈では少し間違ったことを言っていますが、だいたいあってるのであとは調べてください。

もしかしたらこれ単体の記事を書くかもしれません。

実装

ではこれらの概念を使って実際に実装します。

import numpy as np
import cv2

def local_E(spin,i,j,img,img2,img_size):#ローカルな１スピンまわりのエネルギー計算

    E = 0
    for k in range(img_size):
        for l in range(img_size):

            N = img_size**2
            
            #２値(0 or 255)の画像からいちいち直接Jijを計算する。この際0は-1となり、255は1となるように調節
            #予めJijの配列を作成しても良いが、添え字を考えるのがめちゃくちゃ面倒
            #画像が二枚なので手打ちですが、画像自体をリスト化して、for回すのが賢いと思います
            tmp = 2*img[i][j]/255-1 #１枚目画像
            tmp2 = 2*img[k][l]/255-1
            tmp3 = 2*img2[i][j]/255-1 #２枚目画像
            tmp4 = 2*img2[k][l]/255-1
            jij = (tmp*tmp2+tmp3*tmp4)/N #Hebb則からJijを算出
            if i==k and j==l:
                jij = 0.0 #自分自身への結合は０

            E = -1*jij*spin[i][j]*spin[k][l] + E #p番目の画像の全ピクセルに渡ってsummationをとる

    return E

#hyper parameter
img_size = 82 #画像の大きさ(正方形を想定)
border = 105 #２値にする際のボーダーライン、お好みで調整
N = img_size**2 #総ピクセル数
seeds = 11 #乱数シード
sweep = 10 #画像の更新回数(sweep*N回更新)
temp = 0.06 #温度、小さな値だと綺麗に元画像が思い出せるが、計算コストがかかる

#１枚目画像読み込み
img = cv2.imread("pipo.png")
img = cv2.resize(img, (img_size, img_size)) #一応指定したサイズと形にリサイズ
img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) #グレースケール
ret, img = cv2.threshold(img, border, 255, cv2.THRESH_BINARY_INV) #２値化

#二枚目読み込み、上と同様
img2 = cv2.imread("ash.jpg")
img2 = cv2.resize(img2, (img_size, img_size))
img2 = cv2.cvtColor(img2, cv2.COLOR_RGB2GRAY)
ret, img2 = cv2.threshold(img2, border, 255, cv2.THRESH_BINARY_INV)

#２値画像をgazoフォルダに吐き出し(フォルダなどは適当に変えてください)
cv2.imwrite('./gazo/pipo.png', img)
cv2.imwrite('./gazo/ash.png', img2)

#スピン配列(画像を思い出させる用)を作成
spin = np.zeros((img_size,img_size))
tmp_image = np.zeros((img_size,img_size),dtype = int)

#+1-1を適当に配置
for i in range(img_size):
    spin[i] = np.random.choice([-1,1],img_size)

#どのスピン(ピクセル)を更新するかを決めるための配列
xy = np.zeros((img_size),dtype = int)
for i in range(img_size):
    xy[i] = i

np.random.seed(seeds)

for M_i in range(sweep):

    ran = np.random.rand(img_size**2)
    
    k = 0 #画像の吐き出し用

    for i in range(img_size):
        for j in range(img_size):

            #x座標とy座標をランダムにチョイス(普通にやる場合は順番に全部更新する)
            x = np.random.choice(xy,1)[0]
            y = np.random.choice(xy,1)[0]

            E_old = local_E(spin,x,y,img,img2,img_size) #１つのスピンが感じるエネルギーを計算
            spin[x][y] = - spin[x][y] #向きを反転させてみる
            E_new = local_E(spin,x,y,img,img2,img_size) #更新後の１スピン周りのエネルギーを計算
            dE = E_new - E_old #更新前後のエネルギー差を計算
            #本来更新することで系全体のエネルギーが下がるか上がるかを計算しなくてはならない
            #しかし、差分をとるため、注目しているスピンまわりだけ計算すれば十分

            if ran[k] > np.exp(-dE/temp): #乱数が更新採択確率を上回ればreject
                spin[x][y] = - spin[x][y] #元に戻す

            if k%50 == 0: #50回更新するごとに画像を吐き出し(gif用)
                for i in range(img_size):
                    for j in range(img_size):
                        tmp_image[i][j] = int((spin[i][j] * 255 + 255)/2) #スピンを0,255の２値画像に変換

                new_image_path = './gazo/'+str(M_i) +'_'+str(k)+ '.png';
               	cv2.imwrite(new_image_path, tmp_image) #書き込み

            k += 1

ざっとこんなもんです

色々書き足りていないところはありますが、上のプログラムをいい感じに書き直せばもっとたくさんの画像を記憶できるでしょう。

ただ、画像同士がほぼ直交していないとたくさんの画像は記憶できないようです。

結果

上のプログラムだと二枚記憶させているので、実は初期配置によっては

違う人が出てきます。

まとめ

一番簡単な画像記憶をやりましたが、どうでしょうか。

多分説明が足らないところがあると思います。モチベがあればここに付け足すか別の記事で補完したいですね(やらない)

画像記憶以外でもボルツマンマシンですぐできそうなこと見つけたら実装してみようと思います。

藤江杯いってきた話

金沢レジャーランドで開催された第五回藤江杯に行ってまいりました。

詳細はこんな感じ

【第5回ぷよぷよ通大会 藤江杯】
石川県金沢市のレジャーランド藤江本館にて2018年2月24日(土)にぷよぷよ通大会を開催します！
詳細はこちら（https://t.co/zz9H6lVtjU）をご覧ください！#石川ぷよ pic.twitter.com/KXNzgHmoWL

— 藤江ぷよぷよ通対戦会@2/24藤江杯 (@fujieLL_puyo2) 2017年12月24日

地方で開かれていますが、毎回参加人数が50人を超える人気のある大会です。

参加者のほとんどが遠征者という稀有な大会なので、色々な地方のぷよらーをつまみ食いできるいいチャンスです。

僕は今回(も)ぷよらーの方に車を出していただいたので、ワイワイガヤガヤしながら大阪から金沢へーーー

まず、到着後もりもり寿司という「質の高い寿司がたまたま回っているだけ」のお寿司屋さんにいきました。

ついつい3000円使ってしまうものの僕は旅先での予算は常に無限大、財布の紐が緩いどころか底が抜けている状態です。

爽やかな顔で会計を済まし、ACぷよ通のある施設へ出陣します。

最初に初心者大会が行われるわけですが、僕の友人であり犯罪者のぽかりすくんが予選で無双していました。決勝トナメですぐ負けた

僕は出ないことになっていたため一般大会の方に全てをかけていましたが、すんなり負けてしまい真顔で音ゲーコーナーに

しかし、先週普通に落ちて凹んでいた中伝に合格することができ、一気にウキウキモードへ

大会は結局僕の友人であるキウくんが優勝しました。おめでとう

その後は食事会があり、新潟勢を迎え予約していたゆめのゆの大部屋宿泊会場に向かいました。

部屋はこんな感じ

一人一人のスペースも確保されている上に何と言っても「安い」

入館料や部屋代込みで一人当たり3700円程でした。遠征者のみなさん、おすすめです。

「仮眠スペースなら2000円で済むじゃん！」という方々もいらっしゃると思います。

確かに少人数での利用ならその方がいいかも知れませんが、10人集まるならばこのプランがいいと思います。

それはやはり「修学旅行感」に集約されるでしょう。

こんな大部屋に男だけで泊まったらやることは一つ

「恋バナ」

ですよね？？

はい、嘘をつきました。僕はそういったものとは無縁の世界の住人なのでぷよぷよ対戦をしながらツイッターのオタクたちについて話し合っていました。

しかし、お話ししながらゲームしながら広々とスペースを使えるのはアドバンテージですし、テレビもあるので交渉次第ではうまく使えるかも知れません。

ゲーマーの皆さんは是非

翌日は早めに起床し、ジムで大胸筋を鍛えた後、港の食堂で食事をとりました。

これについても色々書きたいですが割愛

最後は忍者寺という全く忍者に関係ない寺に行き、施設内案内ツアーに参加したのですが、キウくんが係員に私語を注意されたところがピークでした。

帰りは車の中でドヒャドヒャ笑いながら運転していたらいつの間にか大阪にいました。

今回の藤江も楽しかったですし、早く就活を終わらせて次の大会にも笑顔で出場したいですね。

では

LSTMっていうオモチャを見付けました

はじめに

これはぼくのインターンメモです。

LSTM

標語的にはリカレントニューラルネットワーク(RNN)の一種でLong short-term memoryの略だとか言われてますね。

この手法の何が他と違うかというと、ある決められた(人の手で入れる)長さのデータを保持しつつ(たまに忘却しつつ)

学習していくというところみたいですね、すごそう。下のリンクの記事は読んだ事ないですが貼ります

qiita.com

何に強いかというと時系列データだとか、文章(前後の繋がりがあるため)でその辺メインで使われてるらしい。

qiita.com
deepage.net

画像認識もいけるみたいだけど、わざわざやる意味はいまいち分かりません。

動かしてみた

株価予想とか腐る程出てくるので他のでやりたい、文章どうこうはあんまり興味がない

んで、今回は

「素数予想」

やります。

CNNで素数判定！みたいな話は結構出てくるんですが、時系列データとして生成するってのはあんまり出てこなかったんで。

まぁ多分やった結果面白くなかったから誰も言ってないないとか、数学的に不可能と証明されているとかそんな感じでしょう。

とりあえず素数をDLしましょう

素数2357: 素数一覧ダウンロード

10000個のやつでやりました、計算重いからね。

毎回思うんですが、こういうcsvファイルってプログラムに読み込ませる以外の用途思いつかないのに、なんで余計なヘッダー付いてるんでしょう。

エクセルで編集はあまりいい思い出が無いので、エディタで消しましょう。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math
from keras.models import Sequential
from keras.models import load_model
from keras.layers import Dense
from keras.layers import LSTM
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error


#データ読み込み
dataframe = pd.read_table('primelist10000.txt')

#変数の抽出ためのデータフレーム作成
df = pd.DataFrame()

#目的変数
df['prime'] = dataframe['Number']#csvのカラム名をぼくが勝手に変えました
max_prime = max(df.prime)#データ範囲での最大値を取得
df.prime = df.prime/max_prime#絶対値0~1でスケール、LSTMはこれが推奨されてるらしい

#説明変数＋上と同じ操作
df['index'] = dataframe.index
max_index = max(df.index)
df.index = df.index/max_index

#ここでnumpy-arrayにする。型は適当に設定
df = df.values
df = df.astype('float32')

#データ分割
train_size = int(len(df) * 0.9)#データの９割で学習、１割でテスト
train, test = df[0:train_size,:], df[train_size:len(df),:]

#LSTMが受け取れる形にしましょう、[データ行数][変数の個数][ルックバックの個数]みたいな形でしか受け取れない
def make_dataset(dataset, look_back=1):
    dataX, dataY = [], []
    for i in range(len(dataset)-look_back-1):
        xset = []
        for j in range(dataset.shape[1]):
            a = dataset[i:(i+look_back), j]
            xset.append(a)
        dataY.append(dataset[i + look_back, 0])
        dataX.append(xset)
    return np.array(dataX), np.array(dataY)


look_back = 10#上で出てきましたが、保持しておくデータ数で今回は適当に１０とおいた
trainX, trainY = make_dataset(train, look_back)#訓練データリシェープ
testX, testY = make_dataset(test, look_back)#テストデータリシェープ

Hidden_Layer = 300#隠れ層の数、多いほどコスト上がるけど多いほどいいとかなんとか
Epochs = 1000#学習回数
Batch_Size = 10#バッチサイズ(まとまった単位で最適化する感じかな？よくわかってなさそう)

#モデル作成
model = Sequential()
model.add(LSTM(Hidden_Layer, input_shape=(testX.shape[1], look_back)))#隠れ層の指定とインプットデータの指定、ルックバックの指定
model.add(Dense(1))#出力の数
model.compile(loss='mean_squared_error', optimizer='adam')#損失関数と最適化法の指定
model.fit(trainX, trainY, epochs=Epochs, batch_size=Batch_Size)#学習開始
model.save("prime_LSTM2.h5")#学習したモデル保存
#model = load_model("prime_LSTM.h5")#読み込む時はこっち

#予測値
trainPredict = model.predict(trainX)
testPredict = model.predict(testX)

#リスケール
trainPredict = trainPredict * max_prime
trainY = trainY * max_prime
testPredict = testPredict * max_prime
testY = testY * max_prime


# 誤差率評価
err_rate_train = (abs(trainPredict-trainY)/trainY)
err_rate_test = (abs(testPredict-testY)/testY)

#以下予測と実際の素数のプロット
test_len = len(testY)
x_pl = np.arange(test_len)
plt.scatter(x_pl,testY,label="actual")
plt.scatter(x_pl,testPredict,label="predict",color="r")
plt.legend()
plt.xlim([0,50])
plt.savefig("ac_pre.png")

#予測された数値を整数として吐き出し
np.savetxt('pre_train.csv',trainPredict.astype("int"),delimiter=',')
np.savetxt('pre_test.csv',(testPredict.astype("int")),delimiter=',')
np.savetxt('test.csv',testY.astype("int"),delimiter=',')

結果

うーん、微妙そうだけど素数が書く軌跡？みたいなのは若干拾えてそう。わからん！

というか予測した値が実際素数かどうか判定しなきゃ意味がないですね！このプロットはなんなんだ！(一応平行移動して心眼で見る事は可能)

そのプログラムは面倒なのでまた今度回しておきます・・・。

まとめ

よくわからん。

というか、int型での推定(離散的な推定)ってそれ用に何かあるのかな

色々改善できそうだけどもうやらなさそうです。

おしまい

Are You Gorilla?

注意　本記事はネタです、ロジックの飛躍や間違いがあります。

お久しぶりです、すしくんです。

前回sym(すわむ)にバーチャルぷよぷよ対戦で勝利してからはや二年————。

未だにリアルぷよぷよ対戦で彼に勝利したことはありません。

ただ現状勝率そのものが変わらない限り人間の寿命では無理、という結論が得られているのである意味仕方ないのかもしれません。

挨拶はこのくらいにして今回のテーマについて話していきたいと思います。

ズバリ

「ぴぽにあゴリラ説の実証」

です。

皆さんも今まで一度は考えた事があると思います。

・・・・





・・・・・・・・・・・





・・・・・・・・・・・・・・・！！！！！！





これでわかって頂けたと思います。

———————————————————————————

しかし、いくら大勢の人が直感的にゴリラだと確信していても全員にアンケートを取るわけにもいきません・・・。

本人に「もしかしてゴリラ？」と聞いたとしても素直に答えてくれるでしょうか？

いずれにせよ客観的に評価する必要があります。




どうしたらいいのか・・・・


最近ハマっているフルグラを貪りながら考えました・・・・


なんでこんなにもフルグラは美味しいのか・・・・


なぜ・・・・・



それは＿＿＿＿











AI・・・・機械学習・・・・！！！！

これがソリューションです。

———————————————————————————

今回主に使用したのはTensorFlowです。

TensorFlow（テンソルフロー）とは、Googleが開発しオープンソースで公開している、機械学習に用いるためのソフトウェアライブラリであるんですよね。

ディープラーニングに対応しており、Googleの各種サービスなどでも広く活用されている。 2017年2月15日に TensorFlow 1.0 がリリースされた^[3][4]。

対応プログラミング言語はC言語、C++、Python、Java、Go^[5]。 対応OSは64ビットのLinux（ただしバイナリ配布はUbuntu用）、macOS、Windows^[6]。ハードウェアは CPU^[7]、NVIDIA GPU^[7]、Google TPU^[8]、Snapdragon Hexagon DSP^[9] などに対応していて、Android Neural Networks API 経由で Android 端末のハードウェアアクセラレータも使用できるんですよね。

今回はpythonで書きました。

・参考記事
TensorFlowでアニメゆるゆりの制作会社を識別する


アルゴリズムはここでは詳しく書きませんが(ちゃんと説明できる気もしない)画像を大量に読み込ませる事で、その人に共通した特徴を学習していくといった感じです。

具体的には２層の畳み込みニューラルネットワークで解析したとか、そういう話もありますがそれは僕のブログの方で書きたいと思います。

そして学習後に、ゴリラ画像を入力しぴぽにあかどうか判定してもらう、という流れで今回は行きたいと思います。(ぴぽにあがゴリラならば、ゴリラを入力した時にぴぽにあと判定される)

また、今回学習に用いたぷよらーは

ぴぽにあ、Ash、キウ、ぽかりす

の顔写真です。

———————————————————————————

まずは画像集めです。これ一番辛い

過去旅行にいった時の画像からサルベージしてきたり、動画から切り抜いてきたり・・・。

集めて切り抜いて・・・・

・・・・・



頭おかしなるで


今回は精度もある程度出したいので、上記4人の顔写真を１００枚程度、ぴぽにあだけ２００枚集めました。

しかしまぁ今回集めたぷよらーは露出が多いので、まだ集めやすかったと思います。（顔は自動で切り抜いていたので、そこまで苦じゃありませんでした。）

そしてこれを学習していきます。



最後のtest accuracyが正答率です。９２パーセントの確率で、ぴぽにあ、Ash、キウ、ぽかりすかどうかを判定できます。そこそこの精度が出ている気がします。(人物判定の相場がわからない)

———————————————————————————

ここからが本番です。

この学習済みモデルを用いて、ぴぽにあが本当にゴリラかどうか判定して行きます。

それでまたゴリラの画像を集めるわけです・・・・




APIの使い方がいまいちよくわからなかったので、手動で集めました。

「このゴリラ・・・・さっき保存した気がする・・・・いや・・・・してない？？」

みたいな謎の苦しみを味わいました。

しかも人間用の顔検出にかけると



このように鼻だけが切り抜かれ、結局手動でやる羽目に・・・。




このゴリラたちをぴぽにあかどうか判定していきます。

———————————————————————————

その結果・・・

集めたゴリラのうち

７５%

が、ぴぽにあであると判定されました。

う〜〜ん、どうなんでしょうこれでは微妙な気がします。

しかしゴリラも色々いるわけで、全てがぴぽにあゴリラであるとは考えられません

そこで「どの程度ぴぽにあなのか」つまり、「ぴぽにあ度」を比較していきます。

今回は

「ぴぽにあと判定された時のぴぽにあである確率の平均」

をぴぽにあ度としました。

まずはぴぽにあ自身のぴぽにあ度を見ましょう。

[ぴぽにあ、Ash、キウ、ぽかりす]

の順番で入力画像が誰なのかの確率が表示されています。




最終的にはぴぽにあ度は一番下の0.997となりました。(理想的なモデルでは１になるべき)


次にゴリラは



0.959となりました。

まとめると

・ぴぽにあのぴぽにあ度は0.997
・ゴリラのぴぽにあ度は0.959

誤差の評価は非常に面倒難しいのでやっていませんが、得られた結論としては

「ぴぽにあはほぼゴリラ」

です。

ゴリラの画像を２００枚くらい集めて同じ計算をすればどんどん1に漸近していくと期待していますが、集めるのが大変なので絶対やりません。誰かゴリラの画像集めて僕に送ってください

———————————————————————————

いかがだったでしょうか？僕は結果にあまり満足していないので、また別のアプローチを考えます。

今回の大きな問題点として「そもそも4人だけの世界で判断している」ので、「4人以外」のラベルを用意し、大量の画像を学習させる必要があります。(ガチのDeep-Learningを用いた高級なアプローチもあるはずです。)

また、おそらくぴぽにあ度は本来ゴリラ全域で平均を取らなければならないとか・・・最初のモデルの精度自体が実は低いとか・・・。

しかしその辺は技術も知識も足りないので、次回かなぁ・・・。









———————————————————————————

オマケ

全ての画像の中で(ぴぽにあの画像も含む)ぴぽにあよりぴぽにあと判断されたゴリラ




Ashと判定されたゴリラ




ぽかりすと判定されたゴリラ




９６パーセントの確率でぴぽにあであると判定された車





キウ君と判定されたゴリラは集めた範囲ではいませんでした・・・。

キウ君は全然ゴリラではないらしい。