## 2012年6月22日(金)
# -*- 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 <- c()
for(iSet in 1:2155){
   activate(iSet);
   hOut <- cbind(hOut, Hidden$out)
}

#中間層の出力を行列の形式で。行方向に15個の中間層の出力、列方向にデータ。

dim(hOut)

par(mfrow=c(1,1))
##相関係数について、おさらい。
##100個の(-1,1)の一様乱数。
x <- runif(n=100,min=-1,max=1) 
hist(x)
y <- x+runif(100,-0.1,0.1) 
plot(x,y)
##相関が高い。
cor(x,y)
title(sprintf("Cor %5.3f", cor(x,y)))
      
y <- runif(n=100,min=-1,max=1) 
plot(x,y)
##相関が低い
cor(x,y)
title(sprintf("Cor %5.3f", cor(x,y)))

y <- -x+0.1*runif(n=100,min=-1) 
plot(x,y)
##相関が高い。
cor(x,y)
title(sprintf("Cor %5.3f", cor(x,y)))

##相関係数の計算cor
cor(hOut[1,],hOut[2,])

##ヒストグラム
hist(hOut)

#中間層の出力のデータに対する、レスポンス。平均と分散。
plot(apply(hOut,2,"mean"))
plot(apply(hOut,2,"var"))

cor(t(hOut))->a
a[row(a)-col(a)>0]
hist(abs(a[row(a)-col(a)>0]),nc=20)

#####「相関係数の辺りの説明用グラフ」書き専用。
dat1 <- read.table("dat1.txt")
plot(dat1,xlim=c(-8.5,2.5),ylim=c(-5,6))

##主成分分析のグラフ用。角度をいれると射影したヒストグラムが出来る。
poi <- function(angle){
  par(mfrow=c(1,2))
  plot(dat1,xlim=c(-8.5,2.5),ylim=c(-5,6))
  abline(a=mean(dat1[,2])-tan(angle*pi/180)*mean(dat1[,1]),b=tan(angle*pi/180)
         ,col="blue")
  d <- matrix(c(cos(angle*pi/180),sin(angle*pi/180)),2,1)##重み付け
  hist( as.matrix(dat1) %*% d, breaks=seq(-10,6,0.5), main="")
  title(paste("angel",angle))
}

for(i in -90:90 ){
  poi(i)
  dev.copy(jpeg, sprintf("D:/tmp/aa%03d.jpg",i+90), width=640, height=480)
  dev.off();
}

sapply(seq(-90,90,1),poi)

par(mfrow=c(1,1))
##主成分分析
par(mfrow=c(1,1))
prcomp(t(hOut))->b

plot((b$sdev)^2,type="b")
abline(h=0)

#累積寄与率
plot(cumsum((b$sdev)^2/sum((b$sdev)^2)),type="b",ylim=c(0,1))

#クラスター解析
dist(hOut)

##2次元平面にプロット
##library(MASS)
##dat.cmd<-cmdscale(as.dist(1.0-cor(t(hOut))))

#15個の中間層の出力の"距離"を総当たりで求める。
hclust(dist(hOut), method="average")
#それに基づいて、階層的クラスタリングを行う。
plot(hclust(dist(hOut), method="average"))

##相関係数の計算cor
cor(hOut[3,],hOut[7,])

##で、相関係数を"距離"と見立てて行うためには、
plot(hclust(as.dist(1.0-cor(t(hOut))), method="average"))


##15次元で分布しているデータを2次元に押しつぶす。
##押し花。
library(MASS)
c <- cmdscale(as.dist(1.0-cor(t(hOut))))
plot(c)


###中間層を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)
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("");
         }
       }
       )
     )