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

##参照ファイルは、
## forPrint7.pdf
## command-7.txt##これ本体。
##・read.table
##・install
##・source
## source("snnAdd.txt")
## source("pat.txt")
## source("snnWeight.txt")

##これをやると、東大の貸出しPCの場合は、
##「作業ディレクトリの変更」をしなくても良い。
setwd("C:\\Users\\student\\Desktop")
##ユーザー>student>デスクトップ
##ニューラルネットの定義
source("snnAdd.txt")
##学習データの設定
source("pat.txt")
##300回の学習結果
source("snnWeight.txt")

##出来上がった、ニューラルネットワークを別の使い方をする。

##エラーの様子。
plot(sapply(1:2155, activate), type="b")
##失敗例について、
abline(h=0.3)
##右辺が比較演算だから、dameは、失敗を示す{True,False}の列。
dame <- sapply(1:2155, activate)>0.3
##{True}に相当する番号のみを選び出す。
(1:2155)[dame]
## 231 389 408 486 835 <<これが失敗した例の通し番号
##それらを、col="red", pch=3で出してみた。
points((1:2155)[dame],
       (sapply(1:2155, activate))[dame],
       col="red", pch=3)

##9塩基の配列は、BufInputのおまけとして付けて有る。
(dimnames(BufInput)[[2]])[dame]
## "CTGGTTCTC" "AGGGGCCAC" "ACAGTGCTG" "CTACAACGG" "CAGGTATAT"

iset <- 230##エラーの少ないやつ。
activate(iset)
Error
(dimnames(BufInput)[[2]])[iset]
(dimnames(BufTarget)[[2]])[iset]


iset <- 231##エラーの多いやつ
activate(iset)
Output$out
## 0.000443492 0.999546532
BufTarget[,iset]
## 1 0
Error
## 0.4995516

##で、この231番の入力配列がどうだったら、エラーが少なくなる。
##つまり、「EI配列らしく」なるのだろうか??

##入力層のデルタを見てみる。
##これは、ほぼlearnと同じ
calDerivative <- function(iset){
  activate(iset);
  Output$delta <<- Output$delta * sapply(Output$out, dsigmoid);
  Hidden$delta <<- drop(t(HiddenOutput$w) %*% Output$delta) * sapply(Hidden$out, dsigmoid);
  ## *sapply(Input$out, dsigmoid);が有ったが、dsigmoidを無視。
  Input$delta <<- drop(t(InputHidden$w) %*% Hidden$delta);
}

calDerivative(iset)

Output$delta
Hidden$delta
matplot(cbind(Input$delta*100,Input$out),type="b")
abline(v=0:9*4+0.5)
##区切りの縦線
text(1:36,0.8,c("A","T","G","C"))##足りないので、循環してくれる
title(sprintf("%d %s", iset, (dimnames(BufInput)[[2]])[iset]))
##凡例
legend("topleft", legend = c("delta","input"),
       col = 1:2, pch=as.character(1:2),
       lty=1, bg = "white")##白抜き

##231番目を書き換える。その前に、保存。
tmpBufInput <- BufInput[,iset]
lines(sapply(Input$delta*10000, sigmoid),type="b",col=3)

##極端に、緑の様にしちゃる。ほぼ、deltaの言いなり。
##sigmoid関数をかけて、[0,1]に規格化している。
BufInput[,iset] <- sapply(Input$delta*10000, sigmoid)
activate(iset)
Error

##やりすぎ??4つの塩基が置換している。
plot(sapply(Input$delta*10000, sigmoid),type="b",col=3)
lines(tmpBufInput,type="b",col=1)
abline(v=0:9*4+0.5)
text(1:9*4-1.5,0.9,1:9,col=2)


##SNPで何とかなるの??
##1~9塩基のうちの1つだけ換えてみる。
mat <- matrix(0,nrow=4,ncol=4)
diag(mat) <- 1 ##対角成分だけに1を
mat

tmp <- tmpBufInput  ##さっき保存した231番目

##2つ目の塩基を、Aに置換する。
##4つの数字を、一度に書き換えている。
tmp[1:4+4] <- as.vector(mat[,1])
tmp  ##変更後
tmpBufInput  ##変更前

ErrorS <- c()
for(ibase in 1:9){
  ##1~9の場所の塩基を、4種類に置換する。
  for(ichange in 1:4){
    ##交換される塩基
    tmp <- tmpBufInput  ##一回貯めて、
    tmp[1:4+(ibase-1)*4] <- as.vector(mat[,ichange])
    BufInput[,iset] <- tmp
    activate(iset)
    ErrorS <- c(ErrorS,Error)
  }
}

plot(ErrorS,type="b")
abline(v=0:9*4+0.5)
lines(tmpBufInput*0.3+0.1, col=2, type="b")

text(1:9*4-1.5,0.2,1:9,col=2)
text(1:36,0.3,c("A","T","G","C"))
##めっけ。6塩基目をAにすれば、良い。

(dimnames(BufInput)[[2]])[iset]
##元々は、"CTGGTTCTC" だった。
##それを、"CTGGTACTC" にすれば、良いってこと。
grep("CTGGTTCTC", (dimnames(BufInput))[[2]])
## 231 元々231番目だった。
grep("CTGGTACTC", (dimnames(BufInput))[[2]])
## 対応する入力データが無い。(ある意味)新規。