#
iz<-sample(seq(1,7291),500)
Nmat<-zip.train[iz,-1]
Nlab<-zip.train[iz,1]
#
center<-function(x) { x<-(x-mean(x))}
#
Nmat<-t(apply(Nmat,1,center))
#
ssp<-svd(Nmat) # SVD on X
ssp2<-svd(t(Nmat)) # SVD on X'
vfr2<-t(Nmat)%*%ssp2$v%*%diag(1/ssp2$d) # getting X'X pcs from XX' pcs
#
plot(ssp$v[,1],vfr2[,1])
abline(0,-1)
plot(ssp$v[,2],vfr2[,2])
abline(0,-1)
# So PCA from X'X can be obtained from PCA on XX'
########
pp<-prcomp(Nmat)
md<-cmdscale(dist(Nmat),2)
plot(md)
plot(Nmat%*%pp$rot[,1],md[,1])
plot(Nmat%*%pp$rot[,2],md[,2])
# PCA and MDS are equivalent if distance is euclidean...
####
####
install.packages("kernlab")
library(kernlab)
#
zeroes<-zip.train[zip.train[,1]==0,-1]
#
ssn<-svd(zeroes)
plot(ssn$d)
plot(ssn$u[,1:2])
pp<-identify(ssn$u[,1:2])
#
par(mfrow=c(2,2))
for (kk in (1:4)) {
image(t(matrix(as.numeric(zeroes[pp[kk],]),16,16,byrow=T))[,16:1])
# highligthing the structure that the 1st components identify and separate
}
ssk<-kpca(zeroes, kernel = "rbfdot", kpar = list(sigma = 0.1),
     features = 50, th = 1e-4)
plot(eig(ssk))
plot(pcv(ssk)[,1:2])
pp<-identify(pcv(ssk)[,1:2])
#
par(mfrow=c(2,2))
for (kk in (1:4)) {
  image(t(matrix(as.numeric(zeroes[pp[kk],]),16,16,byrow=T))[,16:1])
}
####
iz<-sample(seq(1,7291),2000)
Nmat<-zip.train[iz,]
ssn<-svd(Nmat[,-1])
plot(ssn$u[,1:2],pch=as.character(Nmat[,1]),col=gray(Nmat[,1]/10))
#
ssk<-kpca(Nmat[,-1], kernel = "rbfdot", kpar = list(sigma = 0.0005),
          features = 50, th = 1e-4)
plot(pcv(ssk)[,1:2],pch=as.character(Nmat[,1]),col=gray(Nmat[,1]/10))
#
ssk<-kpca(Nmat[,-1], kernel = "polydot", kpar = list(degree = 3),
          features = 50, th = 1e-4)
plot(pcv(ssk)[,1:2],pch=as.character(Nmat[,1]),col=gray(Nmat[,1]/10))
#
ssk<-kpca(Nmat[,-1], kernel = "polydot", kpar = list(degree = 5),
          features = 50, th = 1e-4)
plot(pcv(ssk)[,1:2],pch=as.character(Nmat[,1]),col=gray(Nmat[,1]/10))
#
### tSNE
library(tsne)
ecb = function(x,y){ plot(x,t='n'); text(x,labels=Nmat[,1], col=gray(Nmat[,1]/10)) }
tsne_N = tsne(Nmat[,-1], epoch_callback = ecb, perplexity=25,max_iter=250) # takes a while to run if you run many iterations
# perplexity - bandwidth choice
par(mfrow=c(1,1))
plot(tsne_N,col="white")
text(tsne_N,labels=Nmat[,1],col=gray(Nmat[,1]/10))
# MDS
library(cluster)
library(stats)
ddNbr<-daisy(Nmat[,-1])
ffNbr<-cmdscale(ddNbr,2) # 
plot(ffNbr,pch=as.character(Nmat[,1]),col=gray(Nmat[,1]/10))
#
ddNbr<-as.dist(1-cor(t(Nmat[,-1])))
ffNbr<-cmdscale(ddNbr,2) # 
plot(ffNbr,pch=as.character(Nmat[,1]),col=gray(Nmat[,1]/10))
#
library(vegan)
isoNbr<-isomap(dist(Nmat[,-1]),ndim=2,k=10) # use only 10 nearest neighbors to produce the 2dim representation
plot(isoNbr$points,col=gray(Nmat[,1]/10),pch=as.character(Nmat[,1]))
#
library(lle)
lleNbr<-lle(Nmat[,-1],m=2,k=10,reg=1)
plot(lleNbr$Y[,1:2],col=gray(Nmat[,1]/10),pch=as.character(Nmat[,1]))
#######
install.packages("dimRed") # all the above functions with a wrapper
library(dimRed)
#
lle <- LLE()
emb <- lle@fun(Nmat[,-1], lle@stdpars)
## using embed():
emb2 <- embed(dat, "LLE", knn = 45)
plot(emb, type = "2vars")
plot(emb2, type = "2vars")

######