# -*- mode: R -*-
##EmacsとESS使いの方に。モードラインの設定。
## 2013年7月2日(火曜日)
##
##plot(1:10)
##setwd(choose.dir())
## SVMの判別関数。
x <- seq(-5, 5, by=0.2)
y <- x
g <- function(x1,x2) { 0.2*x1+0.8*x2+0.2 }
z <- outer(x, y, g)
persp(x, y, z, theta = 30, phi = 30)
contour(x, y, z)
abline(h=0,v=0)

##
g<-function(h,l){sqrt((h*log(2*l/h+1)-log(0.001))/l)}
z <- outer(seq(1,100), seq(100,1000,by=100), g)
persp(seq(1,100), seq(100,1000,by=100), z, xlab="h", ylab="#sample", theta = -30, phi = 30)
contour(seq(1,100), seq(100,1000,by=100), z,xlab="h", ylab="#sample")
        

##(学習済みの)SVMの出力(もどき)
## x <- seq(-5, 5, by=0.2)
## y <- x
f <- function(x) {x}
g <- function(x1,x2) { f(0.2*x1+0.8*x2+0.2) }
z <- outer(x, y, g)
persp(x, y, z, theta = 30, phi = 30)
abline(h=0,v=0)
contour(x, y, z)

##ルートの中身は
g<-function(h,l){sqrt((h*log(2*l/h+1)-log(0.001))/l)}
z <- outer(seq(1,100), seq(100,1000,by=100), g)
persp(seq(1,100), seq(100,1000,by=100), z, xlab="h", ylab="#sample", theta = -30, phi = 30)

##|W|を変えてみる
x <- seq(-5, 5, by=0.1)
y <- x
par(mfrow=c(1,2))
i <- 0
for(a in seq(0.1,5,0.1)){
  f <- function(x) { x }
  g <- function(x1,x2) { f(a*(0.2*x1+0.8*x2+0.2)) }
  
  persp(x, y,
        outer(x, y, g),
        theta = 30, phi = 30, zlim=c(-5,5))
  title(sprintf("a=%3.1f",a))

  image(x, y,outer(x, y, g),
        col=terrain.colors(10))
  contour(x, y,outer(x, y, g),
          levels=-100:100, method="simple",add=T)
  dev.copy(jpeg, sprintf("D:/tmp/ch%03d.jpg",i), width=640, height=480)
  dev.off();
  i <- i+1
}

## outerは行列の外積を求める関数。であるが、これはちょっと違う使い方
paste(":","-",")")
paste(":","-",")", sep="")
f <- function(x,y){paste("<",x,", ",y,">", sep="")}
f(3,-5)
x <- 1:5
y <- -2:2
outer(x,y,f)


##hseq <- seq(1.5, 2.3, 0.1) #身長の配列。
##jseq <- seq(0.3, 1.0, 0.1) #ジャンプの配列。30cm〜100cm

hseq <- seq(0, 5.0, 0.1) #身長の配列。
jseq <- seq(0, 5.0, 0.1) #ジャンプの配列。
f <- function(x){ 1/(1+exp(-x))}
g <- function(h,j) { f(1.2*h+0.8*j-2.9) }
z <- outer(hseq, jseq, g) ##身長とジャンプ力の2次元平面の各点に対して、gが適用され、Zが決まる。
contour(hseq, jseq, z)
persp(hseq, jseq, z, theta = 30, phi = 30)  #3Dのグラフを書く

##aが大きくなると、(2,2)の位置での、g()の値が大きくなる。
gga <- function(a) { x1<-2; x2<-2; f(a*(0.2*x1+0.8*x2+0.2)) }
plot(1:5, gga(1:5), ylim=range(c(0, range(gga(1:5)))), type="b")
abline(h=0)

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)
contour(x1,x2,taiseki)
persp(x1,x2,taiseki,theta = 30, phi = 30)


##実際に2次計画問題を解いてみる。
install.packages("e1071", dependencies=TRUE)

##ライブラリの読み込み
library(e1071)
##データの読み込み
source("patnew.dat")
##svmの実行
model <- svm(x, y, kernel="linear")

summary(model)
pred <- predict(model, x)

##サポートベクターに選ばれた入力データは
model$index

##kernel="linear"の場合、平面なので、ニューラルネットワークの1つのセル
##と同じ
apply(model$SV, 2, mean)

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


#予測と本当のラベル
table(pred, y) 

##kernel関数の変更
poi <- function(k){
   model <- svm(x, y, kernel=k, coef0=1);
   summary(model);
   pred <- predict(model, x);
   table(pred, y);
   plot(cmdscale(dist(x)),
     col = as.integer(y),
     pch = c("o","+")[seq(along=y) %in% model$index + 1]);
   title(k)
}
poi("linear")
poi("polynomial") ##coef0=1, dgree=3(default)
poi("radial")

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