# -*- mode: R -*-
##EmacsとESS使いの方に。モードラインの設定。
## 2014年7月18日(金曜日)
# -*- mode: R -*-
##EmacsとESS使いの方に、モードラインの設定。

##plot(sin)
##ディレクトリの変更は、これでもできる。
##setwd(choose.dir())

## EIが677個、Nが1478個。全部で2155個。
source("pat.txt")#これは、dump()コマンドで作成した物、

## #ニューラルネットワークのプログラム
source("snn.txt")

##300回学習済みの重みパラメータをセットする。
source("snnWeight.txt")

##で、中間層の出力は、どうなっているの??
#中間層の出力を行列の形式でhOutに貯める。
hOut <- c()
for(iSet in 1:2155){
  activate(iSet);
   hOut <- cbind(hOut, Hidden$out)
}


##行方向に15個の中間層の出力、列方向にデータ。
dim(hOut)

##15次元で分布しているデータを2次元に押しつぶす。
##押し花。
library(MASS)
c <- cmdscale(as.dist(1.0-cor(t(hOut))))
plot(c,col=2)
text(c[,1]+0.03,c[,2],as.character(1:15))


###中間層を1つにしてみる。
#学習結果
source("snnWeight1.txt")
Hidden
## $net
## [1] 49.89437
## $out
## [1] 1
drop(InputHidden$w)

(1:2155)[sapply(1:2155, activate)>0.1]

##重みの様子を見てみる。
barplot(drop(InputHidden$w),
        col=rep(c("red","green","blue","yellow"),9)
)

dim(InputHidden$w)
matrix(InputHidden$w, nrow=4)
barplot(as.vector(rbind(matrix(InputHidden$w, nrow=4),0)),
        col=c("red","green","blue","yellow","white")
        )


##オマケ、lda()でもやってみよう。
library(MASS)

i <- as.vector(t(BufInput))
im <- matrix(i,nr=2155)
t <- dimnames(BufTarget)[[2]]
b <- lda(im,t)
         
plot(b)
matplot(predict(b,im)$posterior, type="b")
hist(predict(b,im)$x)

b$scaling

barplot(-drop(b$scaling),
   col=rep(c("red","green","blue","yellow"),9)
)


##判別分析でもやってみる
library(MASS)

i <- as.vector(t(BufInput))
im <- matrix(i,nr=2155)
dim(im)
t <- as.vector(t(BufTarget))
tm <- matrix(t,nr=2155)
dim(tm)

b <- lda(im,tm[,2])
plot(b)
##matplotでは、列毎に1本の折れ線グラフを書く。この場合は2本。
matplot(predict(b,im)$posterior)

##hist(predict(b,im)$x)

barplot(drop(b$scaling),
        col=rep(c("red","green","blue","yellow"),9)
)

barplot(as.vector(rbind(matrix(b$scaling, nrow=4),0)),
        col=c("red","green","blue","yellow","white")
        )

##数だけから推定。
##EI境界のもの
apply(im[tm[,1]==1,],2,sum)
##そうじゃないのもの
apply(im[tm[,1]==0,],2,sum)

barplot(apply(im[tm[,1]==1,],2,sum)-apply(im[tm[,1]==0,],2,sum),
        col=rep(c("red","green","blue","yellow"),9))

barplot(as.vector(rbind(-1*matrix(apply(im[tm[,1]==1,],2,sum)-apply(im[tm[,1]==0,],2,sum),
                                  nrow=4),0)),
        col=c("red","green","blue","yellow","white")
        )

##強引に、全ての入力データをクラスタリングしてみる。
##1つめのデータ
im[1,]
## [1] 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0
##2つめのデータ
im[2,]
## [1] 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 1

##差の2乗は、
(im[1,]-im[2,])^2
##その総和、
sum((im[1,]-im[2,])^2)
##総和の最小値は、0.0。これが、ニューラルネットワークでは「2乗誤差」
sum((im[1,]-im[1,])^2)

##つまり、まったく同じ塩基配列だったら、距離は0。2つの塩基配列中に、異
##なる塩基が有れば、距離は増える。このとき、塩基の(違いの)種類によらな
##い。異なる塩基数だけが距離の元。


##全部だと大変
sss <- runif(2155)<0.5
#sss <- runif(2155)<0.2
##plot(hclust(dist(im[sss,])))
hclust(dist(im[sss,]))->hc
##plot(hc)
##ラベルでEIかそうでないかを知る。
##hc$orderで、デンドログラム中の葉が判る。それが==1だと、

plot(hc,
     labels=sapply((tm[sss,1])[hc$order]==1,
       function(x){
         if(x){
           return("*");##EI境界にある
         } else {
           return("");
         }
       }
       )
     )


## SVMの判別関数。
x <- seq(-5, 5, by=0.2)
y <- x
g <- function(x1,x2) { 0.2*x1+0.8*x2+0.2 }
z <- outer(x, y, g)
persp(x, y, z, theta = 30, phi = 30)
contour(x, y, z)
abline(h=0,v=0)

##
g<-function(h,l){sqrt((h*log(2*l/h+1)-log(0.001))/l)}
z <- outer(seq(1,100), seq(100,1000,by=100), g)
persp(seq(1,100), seq(100,1000,by=100), z, xlab="h", ylab="#sample", theta = -30, phi = 30)
contour(seq(1,100), seq(100,1000,by=100), z,xlab="h", ylab="#sample")
        

##(学習済みの)SVMの出力(もどき)
## x <- seq(-5, 5, by=0.2)
## y <- x
f <- function(x) {x}
g <- function(x1,x2) { f(0.2*x1+0.8*x2+0.2) }
z <- outer(x, y, g)
persp(x, y, z, theta = 30, phi = 30)
abline(h=0,v=0)
contour(x, y, z)

##ルートの中身は
g<-function(h,l){sqrt((h*log(2*l/h+1)-log(0.001))/l)}
z <- outer(seq(1,100), seq(100,1000,by=100), g)
persp(seq(1,100), seq(100,1000,by=100), z, xlab="h", ylab="#sample", theta = -30, phi = 30)

##|W|を変えてみる
x <- seq(-5, 5, by=0.1)
y <- x
par(mfrow=c(1,2))
i <- 0
for(a in seq(0.1,5,0.1)){
  f <- function(x) { x }
  g <- function(x1,x2) { f(a*(0.2*x1+0.8*x2+0.2)) }
  
  persp(x, y,
        outer(x, y, g),
        theta = 30, phi = 30, zlim=c(-5,5))
  title(sprintf("a=%3.1f",a))

  image(x, y,outer(x, y, g),
        col=terrain.colors(10))
  contour(x, y,outer(x, y, g),
          levels=-100:100, method="simple",add=T)
  ##dev.copy(jpeg, sprintf("D:/tmp/ch%03d.jpg",i), width=640, height=480);  dev.off();
  i <- i+1
}

## outerは行列の外積を求める関数。であるが、これはちょっと違う使い方
paste(":","-",")")
paste(":","-",")", sep="")
f <- function(x,y){paste("<",x,", ",y,">", sep="")}
f(3,-5)
x <- 1:5
y <- -2:2
outer(x,y,f)


par(mfcol=c(1,2))
##表面積一定の立方体の体積は
x1 <- seq(0,2,by=0.02)
x2 <- x1
## 制約式:x1*x2+x2*x3+x3*x1=3より
## x3*(x1+x2)=3-x1*x2だから
## 体積を求める関数は、x1とx2を用いて
g <- function(x1,x2){x1*x2*(3.0-x1*x2)/(x1+x2)}
taiseki <- outer(x1,x2,g)
contour(x1,x2,taiseki)
persp(x1,x2,taiseki,theta = 30, phi = 30)


##実際に2次計画問題を解いてみる。
install.packages("e1071", dependencies=TRUE)

##ライブラリの読み込み.
##場所として、日本を選ぶ
library(e1071)
##データの読み込み
source("patnew.dat")
##svmの実行
model <- svm(x, y, kernel="linear")

summary(model)
pred <- predict(model, x)

##サポートベクターに選ばれた入力データは
model$index

##kernel="linear"の場合、平面なので、ニューラルネットワークの1つのセル
##と同じ
apply(model$SV, 2, mean)

barplot(apply(-model$SV, 2, mean),
        col=rep(c("red","green","blue","yellow"),9))


#予測と本当のラベル
table(pred, y) 

##kernel関数の変更
poi <- function(k){
   model <- svm(x, y, kernel=k, coef0=1);
   summary(model);
   pred <- predict(model, x);
   table(pred, y);
   plot(cmdscale(dist(x)),
     col = as.integer(y),
     pch = c("o","+")[seq(along=y) %in% model$index + 1]);
   title(k)
}
poi("linear")
poi("polynomial") ##coef0=1, dgree=3(default)
poi("radial")

##cmdscaleは、要素間の違いに基づき、2次元平面に分布させる方法
#distは、要素間の距離行列を計算する。ユークリッド距離。
#col = as.integer(y)で、0(黒)、1(赤)で各点を書く。
#pchで、各点に現れる文字を指定する。
#model$indexは、サポートベクトルを示す
#%in%等を使って、サポートベクトルか否かの区別をしている。match(x, table)