txt = "You say goodbye and I say hello.";txt = tolower(txt)
txt = gsub('[.]', ' .',txt)
words = unlist(strsplit(txt, " "))
uniq.words = unique(words)
uniq.words[1]
which(uniq.words=="say")

n.uniq = length(uniq.words)
n.words = length(words)

# creating corpus
corpus = rep(0,n.words)
corpus = match(words,uniq.words)

## co-occurance matrix
cc = matrix(0,nrow=n.uniq, ncol=n.uniq)
for (i.c in 1:n.uniq){
  loc = which(corpus==i.c)
  L = which(match(uniq.words,words[pmax(loc-1,1)])!='NA')
  R = which(match(uniq.words,words[pmin(loc+1,n.words)])!='NA')
  cc[i.c, c(L,R)]=cc[i.c, c(L,R)]+1
}
diag(cc)=0

## ppmi matrix
ppmi <- function(cc, eps = 1e-8){
  n.uniq = ncol(cc)
  PPMI = matrix(0, n.uniq, n.uniq)
  N = sum(cc)
  r.sum = rowSums(cc)
  pmi = log2(cc*N/outer(r.sum,r.sum))
  PPMI=pmax(0,pmi)
  return(PPMI)
}

# performing SVD
s = svd(PPMI)
plot(s$u[,1],s$u[,2],pch=20,col='red')
text(s$u[,1],s$u[,2],row.names(cc))

# calc. simliarity b/w words
words.sim <- function(word1, word2, eps=1e-8){
  nw1 = word1/(sqrt(sum(word1^2)) + eps)
  nw2 = word2/(sqrt(sum(word2^2)) + eps)
  return(nw1%*%nw2)
}

w1 = cc[which(uniq.words=="you"),]
w2 = cc[which(uniq.words=="i"),]
words.sim(w1,w2)

# reporting most simliar words
most.sim <- function(word, uniq.words, cc, N=5){
  n.uniq = length(uniq.words)
  word2 = cc[which(uniq.words==word),]
  res = data.frame(words = uniq.words, similarity = rep(0,n.uniq))
  top = data.frame(words = rep("",N+1), similarity=rep(0,N+1))
  res$similarity = apply(cc,1, function(x) words.sim(x,word2))
  sort.res = sort(res$similarity, decreasing = T,  index.return=T)
  top$words = res$words[sort.res$ix[1:(N+1)]]
  top$similarity = res$similarity[sort.res$ix[1:(N+1)]]
  self.idx = which(top$words == word)
  return(top[-self.idx,])
} 

## Penn TB
f = file("~/courses/CogMod/CMA2018/ptb.tex")
ptb<-readLines(con=f)
ptb = paste(ptb,"")
words = unlist(strsplit(ptb, " "))
uniq.words = unique(words)
n.uniq = length(uniq.words)
n.words = length(words)

corpus = rep(0,n.words)
corpus = match(words,uniq.words)

cc <- contxt(corpus,uniq.words,words)
PPMI<-ppmi(cc)
library(rsvd)
## 時間がかかります。
s = rsvd(PPMI)

sigmoid.forwd <- function(x){
  return(1/(1+exp(-x)))
}

sigmoid.bckwd <- function(x, dout){
  y = sigmoid.forwd(x)
  return(dout*(1-y)*y)
}

affine.forwd <- function(x, W, b){
  return(x%*%W + matrix(1, nrow = nrow(x), ncol = 1)%*%b)
}

affine.bckwd <- function(x, W, b, dout){
  dx = dout%*%t(W)
  dW = t(x)%*%dout
  db = colSums(dout)
  return(list(dx = dx, dW = dW, db = db))
}

softmax.forwd <- function(x, target){
  max.x = apply(x,1,max)
  C = ncol(x)
  x = x - max.x%*%matrix(1,nrow=1,ncol=C)
  y = exp(x)/rowSums(exp(x))
  delta = 1e-7;
  R = nrow(as.matrix(y))
  return(-sum(target*log(y + delta))/R)
}

softmax.bckwd <- function(x, target,  dout = 1){
  max.x = apply(x, 1, max)
  R = nrow(x)
  C = ncol(x)
  x = x - max.x%*%matrix(1,nrow=1,ncol=C)
  y = exp(x)/rowSums(exp(x))
  return((y-target)/R)
}

softmax.pred <- function(x, target){
  max.x = apply(x,1,max)
  C = ncol(x)
  x = x - max.x%*%matrix(1,nrow=1,ncol=C)
  y = exp(x)/rowSums(exp(x))
  return(y)
}
N = 100; dim = 2; n.class =3;
x = matrix(0,nrow=N*n.class, ncol = dim)
t = matrix(0,nrow=N*n.class, ncol = n.class)

for (i.cl in 1:n.class){
  for (i.N in 1:N){
    radius = i.N/N
    theta = i.cl*4+4*radius+rnorm(1,0,0.2)
    idx = N*(i.cl-1)+i.N
    x[idx,]=c(radius*sin(theta),radius*cos(theta))
    t[idx,i.cl] = 1
  }
}



# spiral data
plot(x[1:100,1],x[1:100,2],pch=20,col='black',cex=2,xlim=c(-1,1),ylim=c(-1,1))
points(x[101:200,1],x[101:200,2],pch=20,col='blue',cex=2)
points(x[201:300,1],x[201:300,2],pch=20,col='green',cex=2)

init.network <- function(n.neurons){
  n.layer = length(n.neurons)
  W = list(); b = list()
  for (i.layer in 1:(n.layer-1)){
    W[[i.layer]] =
      matrix(rnorm(n.neurons[i.layer]*n.neurons[(i.layer+1)],sd = 0.01),nrow=n.neurons[i.layer])
    b[[i.layer]] =  matrix(rnorm(n.neurons[(i.layer+1)],sd = 0.01), nrow = 1)
  }
  return(list(W = W,b = b))
}

train.x = x
train.y = t
params = init.network(c(2,30,30,3))
batch_size = 50; n.iter =50000; lambda =0.3
n.train = nrow(train.x)
loss = rep(0,n.iter)
for (i.iter in 1:n.iter){
  batch_mask = sample(1:n.train, batch_size)
  x.batch = train.x[batch_mask,]
  t.batch = train.y[batch_mask,]
  a1 = affine.forwd(x.batch,params$W[[1]],params$b[[1]])
  z1 = sigmoid.forwd(a1)
  a2 = affine.forwd(z1,params$W[[2]],params$b[[2]])
  z2 = sigmoid.forwd(a2)
  a3 = affine.forwd(z2,params$W[[3]],params$b[[3]])
  z3 = softmax.forwd(a3,t.batch)
  loss[i.iter] = z3
  dwSM = softmax.bckwd(a3, t.batch, 1)
  dwA3 = affine.bckwd(a2,params$W[[3]],params$b[[3]],dwSM)
  dwSG2 = sigmoid.bckwd(a2,dwA3$dx)
  dwA2 = affine.bckwd(a1,params$W[[2]],params$b[[2]],dwSG2)
  dwSG = sigmoid.bckwd(a1,dwA2$dx)
  dwA1 = affine.bckwd(x.batch,params$W[[1]],params$b[[1]],dwSG)
  params$W[[3]] = params$W[[3]] - lambda*dwA3$dW
  params$b[[3]] = params$b[[3]] - lambda*dwA3$db
  params$W[[2]] = params$W[[2]] - lambda*dwA2$dW
  params$b[[2]] = params$b[[2]] - lambda*dwA2$db
  params$W[[1]] = params$W[[1]] - lambda*dwA1$dW
  params$b[[1]] = params$b[[1]] - lambda*dwA1$db
}

# plotting results
plot(loss,type='l', xlab = "trial")

library(plot3D)
M <- mesh(seq(-1,1,length.out = 300),seq(-1,1,length.out = 300))
temp.x = cbind(as.vector(M$x),as.vector(M$y))
a1 = affine.forwd(temp.x,params$W[[1]],params$b[[1]])
z1 = sigmoid.forwd(a1)
a2 = affine.forwd(z1,params$W[[2]],params$b[[2]])
z2 = sigmoid.forwd(a2)
a3 = affine.forwd(z2,params$W[[3]],params$b[[3]])
cl = softmax.pred(a3)
cl.pred = apply(cl,1,which.max)

image(matrix(cl.pred,300,300))
x2= x*0.5+0.5
points(x2[1:100,1],x2[1:100,2],pch=20,col='black',cex=2)
points(x2[101:200,1],x2[101:200,2],pch=20,col='blue',cex=2)
points(x2[201:300,1],x2[201:300,2],pch=20,col='green',cex=2)