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

##h:身長とJ:ジャンプ力を考慮した「バスケットボール選手」認識関数
g <- function(h,j) { 1.2*h+0.8*j-2.9 }
##
hseq <- seq(1.5, 2.3, 0.1) #身長の配列。150cm〜230cm。
jseq <- seq(0.3, 1.0, 0.1) #ジャンプの配列。30cm〜100cm

z <- outer(hseq, jseq, g)
contour(hseq, jseq, z, xlab="height", ylab="jump")
##abline(a, b) 　a、bは切片と傾き。直線 y = a + bx を描きます。 
abline(a=2.9/0.8, b=-1.2/0.8, col=2, lwd=8)



##h:身長とJ:ジャンプ力を考慮した「バスケットボール選手」認識関数
##間違ったバージョン
g <- function(h,j) { 1.0*h-0.5*j-1.45 }
z <- outer(hseq, jseq, g)
contour(hseq, jseq, z, xlab="height", ylab="jump")
##abline(a, b) 　a、bは切片と傾き。直線 y = a + bx を描きます。 
abline(a=-1.45/0.5, b=1.0/0.5, col=2, lwd=8)

g <- function(h,j) { 1.2*h+0.8*j-2.9 }
g <- function(h,j) { 1.0*h-0.5*j-1.45 }


1;
##plot(sin) ##一度、GraphicsDeviceを表示させないと、いけない??
##setwd(choose.dir())
source("snn.txt")

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

##1番目と2番目の入力データは、CAGGTCTGTとAGGGTGAGC
BufInput[,1:2]
##       CAGGTCTGT AGGGTGAGC
##  [1,]         0         1
##  [2,]         0         0
##  [3,]         0         0
##  [4,]         1         0
##  [5,]         1         0
##  [6,]         0         0
##  [7,]         0         1
##  [8,]         0         0
## .............
## [35,]         0         0
## [36,]         0         1

source("snn.txt")
Hidden <- layer(4)
Hidden
Hidden$net <- c(1,3,5,7)
Hidden$net

##1番目と2番目の入力データは、EI境界のデータ
BufTarget[,1:2]
##      EI EI
## [1,]  1  1
## [2,]  0  0

##・数値化のルール
##Aを(1,0,0,0)、Tを(0,1,0,0)、
##Gを(0,0,1,0)、Cを(0,0,0,1)、
##EI境界を(1,0)、無関係な配列を(0,1)

source("snn.txt")
## #ニューラルネットワークのプログラム
## 余分についているのが、
## InputHidden$w <- matrix(runif(length(Input$net)*length(Hidden$net),min=-0.1,max=0.1), length(Hidden$net), length(Input$net));
## Hidden$bias   <- runif(length(Hidden$net),min=-1,max=1);
## など、重みの値を乱数で初期化する。でないと。。。。

length(Input$net)
#[1] 36
length(Hidden$net)
#[1] 15
length(Output$net)
#[1] 2

a <- activate(1)##1番目のデータをフォワードプロパゲーションし、Errorを貯める
a
##[1] 0.07193615
##エラー(2乗誤差)は、0.07193615
Output
##$out
##[1] 0.6053621 0.3633256
##EIに対しては、0.6053621 
##N(EIじゃない)に対しては、0.3633256
target
target-Output$out
(target-Output$out)^2
sum((target-Output$out)^2)
sum((target-Output$out)^2)/2
sum((target-Output$out)^2)/2/2

#2乗誤差が0でない。

e<-c()
for(i in 1:100){
     learn(1);##1番目のデータのみを100回学習する
     e<-c(e,Error);
}

plot(log(e),type="b")

##うーん。完璧!!

hist(sapply(1:100, activate))
#どのデータに対してもErrorが小さい。出来た!!
#本当?
sapply(800:900, activate)
#えっ。何故何故どうして??

plot(sapply(1:2155, activate))
range(sapply(1:2155, activate))

date()
e<-c()
for(i in 1:10){
  for(iSet in order(runif(2155))){
    learn(iSet);
    e<-c(e,Error);
  }
}
date()

##途中の学習している様子。
plot(log(e))

##エラーを行列の形式に変換。各列が、学習回(2155個の学習)のエラー学習回
##でのエラーの平均を求めよう。

m<-matrix(e,nrow=2155)

mean(m[,1])

mean(m[,10])

plot(log(apply(m,2,mean)), typ="b")


##300回の学習結果を貯めて、後で使う。
dump(c("InputHidden","Hidden","HiddenOutput","Output"),"snnWeight.txt")

c("InputHidden","Hidden","HiddenOutput","Output")##が設定される。

##途中の学習している様子。
source("snnWeight.txt")

##エラーの様子。
plot(sapply(1:2155, activate))
hist(sapply(1:2155, activate),nc=300)
hist(log(sapply(1:2155, activate)),nc=300)


##失敗例について、
(1:2155)[sapply(1:2155, activate)>0.1]

(dimnames(BufInput)[[2]])[231]
##[1] "CTGGTTCTC"

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

##で、この時の中間層は、どうなっているの??

hOut <- c()
for(iSet in 1:2155){
   activate(iSet);
   hOut <- cbind(hOut, Hidden$out)
}

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

dim(hOut)

##相関係数の計算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)

#####これは、グラフ書き専用。データはココのDirには無い。
dat1 <- read.table("dat1.txt")
plot(dat1,xlim=c(-8.5,2.5),ylim=c(-5,6))

##判別分析のグラフ用。角度をいれると射影したグラフが出来る。
poi <- function(angle){
  b <- matrix(c(cos(angle*pi/180),sin(angle*pi/180)),2,1)##重み付け
  hist( as.matrix(dat1) %*% b, breaks=seq(-10,6,0.5), main="")
  title(paste("angel",angle))
  system("sleep 1")
}

poi(-90)
poi(-60)
poi(-30)
poi(  0)
poi( 30)
poi( 60)
poi( 90)

                                        #主成分分析
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"))



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

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


##オマケ、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))