# -*- mode: R -*-
##EmacsとESS使いの方に。モードラインの設定。
## 2016年10月21日 (金)


##文字イメージデータの読み込み
mojiDat <- as.matrix(read.table("t10k-images.txt"))/255
range(mojiDat)##0-255
dim(mojiDat)
## [1] 10000   784 ;;10000個の文字データ、28x28のます。784=28*28

mojiLab <- unlist(read.table("t10k-labels.txt"))
length(mojiLab)
## [1] 10000 ;;10000個の文字データのラベル。0~9
table(mojiLab)
##   0    1    2    3    4    5    6    7    8    9 
## 980 1135 1032 1010  982  892  958 1028  974 1009 

for( i in 1:100 ){
  ##i番目の文字データ(728の実数)を読み込んで。
  tmp <- matrix(mojiDat[i,],ncol=28)
  ##縦方向が下から積みあがっているので、ちょっと換える。
  image(x=1:28, y=1:28, z=tmp[,28:1],col=gray(seq(0,1,len=256)))
  par(new=T)
  contour(x=1:28, y=1:28, z=tmp[,28:1],col=2)
  title(sprintf("%d  (%d bannme)", mojiLab[i], i))
}


source("snnMoji.txt")
##実は前回まで使っていたのと、全く同じのニューラルネット
##違いはニューロンの数の定義のところのみ。
## # network definition;
## # layers
## Input <- layer(28*28);##入力層は、28*28
## Hidden <- layer(30);##中間層は、30個。
## Output <- layer(10);##出力層は、0~9を当てるので、10個。

dim(mojiDat)
## [1] 10000   784 ;;10000個の文字データ、28x28のます。784=28*28

##入力バッファの形式は、
##BufInput[<入力層の個数分>,<サンプルの個数分>]。なので、転置。
BufInput <- t(mojiDat)
dim(BufInput)


##出力バッファの形式は、
##BufTarget[<出力層の個数分>,<サンプルの個数分>]。
BufTarget <- matrix(0, nrow=10,ncol=dim(BufInput)[2])
dim(BufTarget)##全部、0にしたっつ。
BufTarget[,1:4]
mojiLab[1:4]

for( i in seq(dim(BufInput)[2]) ){
  if(mojiLab[i]==0){##おバカ。「0」の時に10にする
    BufTarget[10,i] <- 1; 
  } else {
    BufTarget[mojiLab[i],i] <- 1
  }
}
BufTarget[,1:4]
mojiLab[1:4]

##以上で準備完了っ。

e<-c()
gcc(0)
##いざ学習。
date()
for(i in 1:10){
  for(iSet in order(seq(dim(BufTarget)[2]))){
    learn(iSet);
    if(iSet%%100==0){print(iSet)};
    e<-c(e,Error);
  }
}
date()

aa <- seq(dim(BufInput)[2])
aa[sapply(aa, activate)>0.1]
plot(sapply(aa, activate))

i <- 592
activate(i)
matplot(cbind(Output$out,target), type="b")

##i番目の文字データ(728の実数)を読み込んで。
tmp <- matrix(mojiDat[i,],ncol=28)
##縦方向が下から積みあがっているので、ちょっと換える。
image(x=1:28, y=1:28, z=tmp[,28:1],col=gray(seq(0,1,len=256)))
par(new=T)
contour(x=1:28, y=1:28, z=tmp[,28:1],col=2)
title(sprintf("%d  (%d bannme)", mojiLab[i], i))



##ちょっとズレ、「1」を例にして。
par(mfrow=c(1,1))
for( i in seq(dim(BufInput)[2])[1:1000] ){
  if(mojiLab[i]==1){
    ##i番目の文字データ(728の実数)を読み込んで。
    tmp <- matrix(mojiDat[i,],ncol=28)
    ##縦方向が下から積みあがっているので、ちょっと換える。
    image(x=1:28, y=1:28, z=tmp[,28:1],col=gray(seq(0,1,len=256)))
    par(new=T)
    contour(x=1:28, y=1:28, z=tmp[,28:1],col=2)
    title(sprintf("%d  (%d bannme)", mojiLab[i], i))
    ##画像のモーメントを計算する。
    m00 <- sum(tmp)
    bb1 <- matrix(1:28,ncol=28,nrow=28)
    bb1t <- t(bb1)

    ##生の1次モーメントを計算して
    m10 <- sum(tmp*bb1)
    m01 <- sum(tmp*bb1t)

    ##生の2次モーメントを計算して
    bb2 <- matrix((1:28)^2,ncol=28,nrow=28)
    m20 <- sum(tmp*bb2)
    bb2t <- t(bb2)
    m02 <- sum(tmp*bb2t)
    m11 <- sum(tmp*bb1*bb1t)

    ##重心周りの2次モーメントに変換して
    u20 <- m20/m00 - (m10/m00)^2
    u02 <- m02/m00 - (m01/m00)^2
    u11 <- m11/m00 - (m10/m00)*(m01/m00)

    b=1/(1/2*atan(2*u11/(u20-u02)))##傾きを計算して、
    y0 <- m01/m00
    x0 <- m10/m00

    abline(v=m10/m00, col="white")
    abline(h=m01/m00, col="white")
    ##(x0,y0)を通り、傾きｂの直線。
    ##ablineは、「a+b*x」である。
    abline(a=y0-b*x0, b=1/(1/2*atan(2*u11/(u20-u02))),col="yellow",lwd=3)
  }
}