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

##4本の直線に囲い込まれた領域
x <- seq(-30, 30) ##例によって、-30〜+30の配列を作る
y <- x
f <- function(x){ 1/(1+exp(-x))}##Sigmoid関数

par(mfrow=c(1,1))
plot(x,f(x),type="b")

g1 <- function(x1,x2) { f(0.8*x1-0.5*x2+10) }
g2 <- function(x1,x2) { f(x2+10) }
g3 <- function(x1,x2) { f(-x1+10) }
g4 <- function(x1,x2) { f(-x2+10) }
g <- function(x1,x2) { f(5*(g1(x1,x2)+g2(x1,x2)+g3(x1,x2)+g4(x1,x2)-3.5))}
                            
z <- outer(x, y, g)

par(mfrow=c(1,2))
persp(x, y, z, theta = 30, phi = 30)

##contour(x, y, z)
image(x, y, z, col = terrain.colors(10))
contour(x, y, z, add = TRUE, col = "black")

abline(a=20,b=1.6,col=5,lwd=3) ## y=a+b*x
abline(v=+10,col=4,lwd=3)
abline(h=-10,col=2,lwd=3)
abline(h=+10,col=3,lwd=3)

##################
##色々とパラメータを変えてみたGifアニメ
##保存先のDIRを決める。
setwd("D:/tmp")
x <- seq(-30, 30) ##例によって、-30〜+30の配列を作る
y <- x
f <- function(x){ 1/(1+exp(-x))}##Sigmoid関数

par(mfrow=c(1,1))
plot(x,f(x),type="b")

g1 <- function(x1,x2) { f(0.8*x1-0.5*x2+10) }
g2 <- function(x1,x2) { f(x2+10) }
g3 <- function(x1,x2) { f(-x1+10) }
g4 <- function(x1,x2) { f(-x2+10) }

par(mfrow=c(1,2))
i <- 0
for( a in c(seq(0.1,7,by=0.1),seq(6.9,0.2,by=-0.1))){
  ##g <- function(x1,x2) { f(5*(g1(x1,x2)+g2(x1,x2)+g3(x1,x2)+g4(x1,x2)-3.5))}
  g <- function(x1,x2) { f(a*(g1(x1,x2)+g2(x1,x2)+g3(x1,x2)+g4(x1,x2)-3.5))}
  z <- outer(x, y, g)
  persp(x, y, z, theta = 30, phi = 30)
  image(x, y, z, col = terrain.colors(10))
  contour(x, y, z, add = TRUE, col = "black")
  title(sprintf("%5.3f",a))
  abline(a=20,b=1.6,col=5,lwd=3) ## y=a+b*x
  abline(v=+10,col=4,lwd=3)
  abline(h=-10,col=2,lwd=3)
  abline(h=+10,col=3,lwd=3)
  ##ファイル名を決めて書きだす。
  dev.copy(jpeg, sprintf("D:/tmp/ch%03d.jpg",i), width=640, height=480)
  dev.off();
  i <- i+1
}
1;

##ココからは、Windowsの「ImageMagick」をつかってアニメーションに変換
##pc[d:tmp]% cd D:/tmp
##pc[d:tmp]% convert -delay 10 -loop 0 ch*.jpg animation.gif
##convertは、ImageMagickに含まれるコマンド
##-delay 10::各コマ間のdelay。単位はミリ秒。
##-loop 0:: ループ回数の指定。「0」は、無限ループ。
##最終的に、d:/tmp/animation.gifが出来る。

plot(1:10)
########
x <- seq(-30, 30) ##例によって、-30〜+30の配列を作る
y <- x

##ちょっと細工した。
ff <- function(x,a){ 1/a*(1/(1+exp(-x*a))-1/2)}##Sigmoid関数

par(mfrow=c(1,1))
bb <- 0
for( a in seq(0.1,10,by=0.01)){
  bb <- cbind(bb,ff(x,a))
}
matplot(x, bb, type="l")
lines(x,ff(x,1),type="l",lw=3,col="green")

gg1 <- function(x1,x2) { ff(0.8*x1-0.5*x2+10,a) }
gg2 <- function(x1,x2) { ff(x2+10,a) }
gg3 <- function(x1,x2) { ff(-x1+10,a) }
gg4 <- function(x1,x2) { ff(-x2+10,a) }

a <- 0.1
a <- 0.05

g <- function(x1,x2) { f(5*(gg1(x1,x2)+gg2(x1,x2)+gg3(x1,x2)+gg4(x1,x2)-3.5))}
i <- 0
for( a in seq(0.01,0.2,by=0.002)){
  z <- outer(x, y, g)
  par(mfrow=c(1,2))
  persp(x, y, z, theta = 30, phi = 30)
  
  ##contour(x, y, z)
  image(x, y, z, col = terrain.colors(10))
  contour(x, y, z, add = TRUE, col = "black")
  
  abline(a=20,b=1.6,col=5,lwd=3) ## y=a+b*x
  abline(v=+10,col=4,lwd=3)
  abline(h=-10,col=2,lwd=3)
  abline(h=+10,col=3,lwd=3)

  dev.copy(jpeg, sprintf("D:/tmp/zz%03d.jpg",i), width=640, height=480)
  dev.off();
  i <- i+1
}


###########
##もしも、シグモイド関数じゃなかったら。。。
##

f <- function(x){x}##単なるリニア関数

par(mfrow=c(1,1))
plot(x,f(x),type="b")

z <- outer(x, y, g)
par(mfrow=c(1,2))
persp(x, y, z, theta = 30, phi = 30)

##contour(x, y, z)
image(x, y, z, col = terrain.colors(10))
contour(x, y, z, add = TRUE, col = "black")

abline(a=20,b=1.6,col=5,lwd=3) ## y=a+b*x
abline(v=+10,col=4,lwd=3)
abline(h=-10,col=2,lwd=3)
abline(h=+10,col=3,lwd=3)


########
f <- function(x){ifelse(x>0,1,0)}##閾値関数

par(mfrow=c(1,1))
plot(x,f(x),type="b")

z <- outer(x, y, g)
par(mfrow=c(1,2))
persp(x, y, z, theta = 30, phi = 30)

##contour(x, y, z)
image(x, y, z, col = terrain.colors(10))
contour(x, y, z, add = TRUE, col = "black")

abline(a=20,b=1.6,col=5,lwd=3) ## y=a+b*x
abline(v=+10,col=4,lwd=3)
abline(h=-10,col=2,lwd=3)
abline(h=+10,col=3,lwd=3)

######
f <- function(x){ifelse(x>10, 1, ifelse(x>-10, 0.05*x+0.5, 0))}##区分線形関数

par(mfrow=c(1,1))
plot(x,f(x),type="b")

z <- outer(x, y, g)
par(mfrow=c(1,2))
persp(x, y, z, theta = 30, phi = 30)

##contour(x, y, z)
image(x, y, z, col = terrain.colors(10))
contour(x, y, z, add = TRUE, col = "black")

abline(a=20,b=1.6,col=5,lwd=3) ## y=a+b*x
abline(v=+10,col=4,lwd=3)
abline(h=-10,col=2,lwd=3)
abline(h=+10,col=3,lwd=3)

#################
##ホームページから、たどっていって、data1.txtとdata2.txtをDownLoadする。
##読み込み先ディレクトリの指定。[File]->[Change dir]で変更する。
##choose.dir()
##
choose.dir()
setwd(choose.dir())
##getwd()

dat1 <- read.table("dat1.txt")
dat2 <- read.table("dat2.txt")

##2つのデータをファイルからテーブル形式で読み取る。
##各々1000個のデータからなっている。

dim(dat1)  #行列形式。次元を調べる。
##[1] 1000    2
##これは、1000行*2列。各行は、(x座標,y座標)の形式。

par(mfrow=c(1,1))
plot(dat1) #とやっても出来るが、
plot(dat1,xlim=c(-10,10),ylim=c(-10,10),col="red") #でも出来る
##xlimは、x軸の最大値と最小値の指定。
##ylimは、y軸の最大値と最小値の指定。
##colは色の指定。
##dat2と重ねて見たい。色は青だ。
plot(dat2,xlim=c(-10,10),ylim=c(-10,10),col="blue")
##とすると、前のプロットが消えてしまった。
##作図の関数は、書く前に画面をきれいにしてから、描画する。
##で、重ねて書くためには、作図関数を呼び出す前に

#「もう、画面はきれいにクリアされている」(画面をクリアしないでね)、とオ
#マジナイをかけてから描画する。
#そのための関数がpar(new=T)

plot(dat1,xlim=c(-10,10),ylim=c(-10,10),col="red")
par(new=T)
plot(dat2,xlim=c(-10,10),ylim=c(-10,10),col="blue")

##もう一つのやり方。
plot(dat1,xlim=c(-10,10),ylim=c(-10,10),col="red")
##点だけを書き足す。昔の座標軸が保存されている。
points(dat2,col="blue")


###memo******
##判別分析のグラフ用。以下の2つはすごいナイーブ
a <- hist((rbind(dat1,dat2))[,1], nc=30, plot=F)
hist(dat1[,1],breaks=a$breaks,col="red")
par(new=T)
hist(dat2[,1],breaks=a$breaks,col="blue")

a <- hist((rbind(dat1,dat2))[,2], nc=30, plot=F)
hist(dat1[,2],breaks=a$breaks,col="red")
par(new=T)
hist(dat2[,2],breaks=a$breaks,col="blue")

##判別分析のグラフ用。角度をいれると射影したグラフが出来る。
poi <- function(angle){
  angle
  b <- matrix(c(cos(angle*pi/180),sin(angle*pi/180)),2,1)##重み付け
  a <- hist((as.matrix(rbind(dat1,dat2)))%*%b, nc=30, plot=F)
  ##countの最大値が知りたいだけ。
  c <- hist(as.matrix(dat1)%*%b,breaks=a$breaks,plot=F)
  yylim <- c(0,range(c$counts)[2]*1.1)
  hist(as.matrix(dat1)%*%b,breaks=a$breaks,col="red",main="",ylim=yylim)
  par(new=T)
  hist(as.matrix(dat2)%*%b,breaks=a$breaks,col="blue",main="",ylim=yylim)
  title(paste("angel",angle))
  system("sleep 1")
}
poi(0)
poi(-30)
poi(-60)

sapply(seq(0,360,10),poi)

help(hist)

for(angle in seq(0,360,2)){
  poi(angle)
  dev.copy(jpeg, sprintf("D:/tmp/ldach%03d.jpg",angle), width=640, height=480)
  dev.off();
}


lda()
library(MASS)
help(lda)

##("red","red","red","red","red","red","red",
## "blue","blue","blue","blue","blue","blue","blue")
##が欲しい。

rep("red",7)
rep("blue",7)
rep("red",7),rep("blue",7)

c(rep("red",7),rep("blue",7))
rep(c("red","blue"),3)
rep(c("red","blue"),c(3,5))
rep(c("red","blue"),c(7,7))
rep(c("red","blue"),each=7) 

lda(rbind(dat1,dat2),rep(c("red","blue"),each=1000))

plot(dat1,xlim=c(-10,10),ylim=c(-10,10),col="red")
points(dat2,col="blue")
abline(0,-0.7467589/-1.4874124)


#Coefficients of linear discriminants:
#          LD1
#V1 -1.4874124
#V2 -0.7467589
#help(lda)によると、上記の値は、scalingって所に入っている。
#呼び出すためには、
lda(rbind(dat1,dat2),rep(c("red","blue"),each=1000))->b
plot(b)

(b$scaling)[1]

#分割線の傾きは、(b$scaling)[2]/(b$scaling)[1]
abline(0, (b$scaling)[2]/(b$scaling)[1],col="green")

## dat1からの10個のデータを抜き出す、
## dat2からの10個のデータを抜き出す、
## 抜き出した2つを合体する

## dat1の1行一列目だけ
dat1[1,1]
## dat1の1行目だけ。指定しないと全部の列。
dat1[1,]
## dat1の1列目だけ。指定しないと全部の行。
dat1[,1]
## dat1から1〜10行目だけ。
dat1[c(1,2,3,4,5,6,7,8,9,10),]
dat1[1:10,]

## dat2からの10個のデータを抜き出す、
dat2[1:10,]

## 抜き出した2つを合体する
rbind(dat1[1:10,],dat2[1:10,])

predict(b,rbind(dat1[1:10,],dat2[1:10,]))

predict(b,rbind(dat1[1:10,],dat2[1:10,]))$posterior

hist(predict(b,rbind(dat1,dat2))$x)


##もうちょっと真面目に
plot(dat1+5*runif(2000),xlim=c(-10,10),ylim=c(-10,10),col="red")
points(dat2+5*runif(2000),xlim=c(-10,10),ylim=c(-10,10),col="blue")

hist(predict(b,rbind(dat1+5*runif(2000),dat2+5*runif(2000)))$x, nc=20)->c
hist(predict(b,dat1+5*runif(2000))$x, breaks=c$breaks, col="red");
par(new=T)
hist(predict(b,dat2+5*runif(2000))$x, breaks=c$breaks, col="blue");

##ホームページから、"snn.txt"をDownloadして、Rに貼り付ける。と、ニュー
##ラルネットワークの定義が出来る。

##もしくは読み込む。
source("snn.txt")

##"snn.txt"には、色々な設定が有る。
Eta
Alpha
Maxv
Minv

#sigmoid関数の定義。 
#a<-seq(1:3)と似ている。 aという変数に、(1,2,3)というベクトルを割り当てる。
# sigmoidという変数に、{…}で、定義した関数(function)を割り当てる。
#return()関数は、関数の終了と戻り値の指定。
## sigmoid <- function( net ){#;
##     return((Maxv-Minv)/(1+exp(-net)) + Minv);
## }
sigmoid

sigmoid(4)
plot(seq(-5,5,0.5),sigmoid(seq(-5,5,0.5)),type="b")
abline(v=0)
abline(h=0)

Input
Hidden
Output

matrix(1:9,3,3)%*%c(1,2,3)
drop(matrix(1:9,3,3)%*%c(1,2,3) ) 


#########