## 2017年10月 6日(金曜日)
# -*- mode: R -*-
##EmacsとESS使いの方に、モードラインの設定。
##Up忘れて申し訳ございませんでした。　<(_ _)>

## ##参照ファイルは、
## forPrint*.pdf
## command.txt
##・read.table
##・source
## source("pat.txt")#これは、dump()コマンドで作成した物、
## source("snn.txt")
## source("snnWeight.txt")
## source("patnew.dat")
##・install
## install.packages("e1071_1.6-7.zip", repos=NULL)
##

##plot(sin)
##setwd(choose.dir())

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

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

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

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

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

(dimnames(BufInput)[[2]])[sapply(1:2155, activate)>0.1]
## "CTGGTTCTC"
## "AGGGGCCAC"
## "ACAGTGCTG"
## "CTACAACGG"
## "AAGGTGCCA"
## "CAGGTATAT"
## "AAGGTGCCA"

BufTarget[,sapply(1:2155, activate)>0.1]
##      EI EI EI EI EI N N
## [1,]  1  1  1  1  1 0 0
## [2,]  0  0  0  0  0 1 1

##593はEI,1944はN。でも同じ配列
(dimnames(BufInput)[[2]])[593]
(dimnames(BufInput)[[2]])[1944]

##実際に使うには、SlidingWindow
##1塩基の変換(塩基を実数に)
zz<-function(x){switch(EXPR = x,
                       "A"=c(1,0,0,0),
                       "T"=c(0,1,0,0),
                       "G"=c(0,0,1,0),
                       "C"=c(0,0,0,1),
                       stop("hen"))}##これは、どれにも当てはまらないときに評価される。
zz("A")
zz("T")
zz("N")


##部分配列の切り出しには、substring
##substring("元の文字列",開始位置,終了位置)
substring("ATGC",1,1)
##{開始を1、終了位置を1}とすると、"A"が取れる。
substring("ATGCATGG",2,5)
##{開始を2、終了位置を5}とすると、
#[1] "TGCA"
substring("ATGCATGG",1:3,5)
##{開始を1:3、終了位置をc(5,5,5)}としたのと同じ。
##長さをそろえてくれた。
#[1] "ATGCA" "TGCA"  "GCA" 

substring("ATCC",1:4,1:4)
##{開始を1、終了位置を1}、
##{開始を2、終了位置を2}、
##{開始を3、終了位置を3}、
##{開始を4、終了位置を4}
##[1] "A" "T" "C" "C"

##長い配列を1つ1つに分解できた。ので、
##各々にzz関数に適用(apply)する。
sapply(substring("ATGC",1:4,1:4),zz)
##      A T C C
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 0 0
## [4,] 0 0 1 1

#で欲しいのは、36個の{0,1}のながーいベクトルへの型なので、
#行列になってしまったのを、ベクトルに(無理やり)変換
as.vector(sapply(substring("ATCC",1:4,1:4),zz))
## [1] 1 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1

##ひっくるめて関数化。
##nchar(x)は、文字列の長さ
yy<-function(str){
  as.vector(sapply(substring(str,1:nchar(str),1:nchar(str)),zz))};
##yy("TTCAGCGGCCTCAGCCTGCCTGTCTCC")
yy("TTCAGCGG")


##長いテスト文字列が有った。
##EI,    ATRINS-DONOR-521,
##CCAGCTGCATCACAGGAGGCCAGCGAGCAGGTCTGTTCCAAGGGCCTTCGAGCCAGTCTG
testSeq <- "CCAGCTGCATCACAGGAGGCCAGCGAGCAGGTCTGTTCCAAGGGCCTTCGAGCCAGTCTG"
nchar(testSeq)## 全部で60文字
##9文字ずつに切って、{0,1}数字列に変換
sapply(substring(testSeq,1:(nchar(testSeq)-8), 9:nchar(testSeq)), yy)
1:(nchar(testSeq)-8)
9:nchar(testSeq)


rm(BufInput)##一度、クリアにする。
##BufTargetは、学習しないので設定しない。
BufInput <-
  sapply(substring(testSeq,1:(nchar(testSeq)-8), 9:nchar(testSeq)), yy)

resp <- c()##出力結果を保持する
for(iSet in 1:(nchar(testSeq)-8)){
  activate(iSet)
  Output$out
  resp <- rbind(resp, Output$out)
}
matplot(resp, type="b")

##強引に、全ての入力データをクラスタリングしてみる。
source("pat.txt")#これは、dump()コマンドで作成した物、
BufInput[,1]


##1つめのデータ
BufInput[,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つめのデータ
BufInput[,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乗は、
(BufInput[,1]-BufInput[,2])^2
##その総和、
sum((BufInput[,1]-BufInput[,2])^2)
##総和の最小値は、0.0。これが、ニューラルネットワークでは「2乗誤差」
sum((BufInput[,1]-BufInput[,1])^2)

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


dim(BufInput)
##全部だと大変
sss <- runif(2155)<0.5
sum(sss+0)
##runifはrとunif
##rは乱数を意味する。unifは[0.0, 1.0]の一様分布
head(sss)
hist(sss+0)
hc <- hclust(dist(t(BufInput[,sss])))
plot(hc)

##plot(hc)
##ラベルでEIかそうでないかを知る。
##hc$orderで、デンドログラム中の葉の順番が判る。
dimnames(BufTarget)[[2]]##で、おまけに付けたEIかNかが判る。
((dimnames(BufTarget))[[2]])[sss]##上記で作った、選択基準を適用。sssはTもしくはF
(((dimnames(BufTarget))[[2]])[sss])[hc$order]##クラスタリングで並べられた順。
(((dimnames(BufTarget))[[2]])[sss])[hc$order]=="EI"##ラベルが、EIか??
ifelse((((dimnames(BufTarget))[[2]])[sss])[hc$order]=="EI","*","")##そうだったら*を。
plot(hc,
     labels=ifelse((((dimnames(BufTarget))[[2]])[sss])[hc$order]=="EI","*",""))



##こっからSVM
##こっからSVM
##こっからSVM

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)
image(x1,x2,taiseki,col=terrain.colors(30))
par(new=T)
contour(x1,x2,taiseki)
persp(x1,x2,taiseki,theta = 30, phi = 30)

##実際に2次計画問題を解いてみる。
##setwd(choose.dir())
install.packages("e1071_1.6-7.zip", repos=NULL)
##install.packages("e1071", dependencies=TRUE)
library(e1071)

##データの読み込み
source("patnew.dat")
seqS <- rownames(x)##配列を保存する。
rownames(x) <- c()##rowNameを空にする。
y <- ifelse(y==1,"EI","N")##ラベルの付け替えをする。
y <- as.factor(y)##ファクタに変換する。

##svmの実行
model <- svm(x, y, kernel="linear")
summary(model)

##予測する。
pred <- predict(model, x)
##サポートベクターに選ばれた入力データはkernel="linear"の場合、平面なの
##で、ニューラルネットワークの1つのセルと同じ
table(pred, y)

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

##予測と本当のラベルの予測率。
(aaa <- table(pred, y))
##対角要素が正解。
sum(diag(aaa))/sum(aaa)

##kernel関数の変更
poi <- function(k){
   model <- svm(x, y, kernel=k, coef0=1);
   summary(model);
   pred <- predict(model, x);
   aaa <- table(pred, y);
   plot(cmdscale(dist(x)),
     col = as.integer(y),##Labelで色を変える。
     pch = c("o","+")[seq(along=y) %in% model$index + 1]);
   title(sprintf("%s correct %2.2g%%", k, sum(diag(aaa))/sum(aaa)*100))
}

poi("linear")
poi("polynomial") ##coef0=1, dgree=3(default)
poi("radial")

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


##マージンを最大化しておくと将来役立つ。
##これが、「未学習データに対しても良い」の理屈。
maru <- matrix(
               c(13, 72,
                 11, 35,
                 37, 57,
                 23, 41,
                 26, 21,
                 40, 68,
                 43, 50,
                 68, 75,
                 89, 77),
               ncol=2 ,byrow=T)
batsu <- matrix(
                c(68, 50,
                  77, 43,
                  88, 46,
                  62, 27,
                  78, 29,
                  62, 31,
                  55, 13,
                  89, 16),
                ncol=2 ,byrow=T)

##svmのマージン最大化の説明
plot(rbind(maru,batsu),type="n")
##pch=21は色ありの●。
points(maru[,],pch=21,col=1,cex=2,lw=3,bg="orange")
points(batsu[,],pch=4,col=2,cex=2,lw=3)

col2rgb("orange")/255
plot(rbind(maru,batsu),type="n")
##21は色ありの●。
for( i in 1:1000 ){
  points(maru[,]+5*rnorm(length(as.numeric(maru))),
         col=NULL, bg=rgb(1, 0.64, 0, alpha=0.1), pch=21, cex=2)
}
points(maru[,],pch=21,col=1,cex=2,bg="orange",lw=3)

for( i in 1:1000 ){
  points(batsu[,]+5*rnorm(length(as.numeric(batsu))),
         pch=4,col=1,cex=2)
}
points(batsu[,],pch=4,col=2,cex=2,lw=2)