library("NMF")
# a subset of handwritten digits
N<-500
# sample 500 digits to try NMF on
it<-sample(seq(1,7291),500)
K<-10
aa<-nmf(Nmat[it,]-min(Nmat),K)
coefmap(aa)
# this is what the codebook looks like (each is a raster scan of a 16*16 image)
plot(basis(fit(aa))[,3:4],col=as.numeric(Numbers[it,1])+1,
             ,pch=as.character(Numbers[it,1]))
# plotting basis functions 1:2 - separates which digits?
#
# Try a few of the digits and look what the K-rank rep captures
image(matrix(Nmat[it[2],],16,16,byrow=T))
image(matrix(fitted(aa)[2,],16,16,byrow=T))
#####
# Try with different K
#
summary(aa)
# percent of variance explained
1-sum(((Nmat[it,]-min(Nmat))-fitted(aa))^2)/sum((Nmat[it,]-min(Nmat))^2)
#
K<-3
aa<-nmf(Nmat[it,]-min(Nmat),K)
coefmap(aa)
summary(aa)
1-sum(((Nmat[it,]-min(Nmat))-fitted(aa))^2)/sum((Nmat[it,]-min(Nmat))^2)
#
### different algorithms
K<-3
aa<-nmf(Nmat[it,]-min(Nmat),K,"lee")
coefmap(aa)
summary(aa)
1-sum(((Nmat[it,]-min(Nmat))-fitted(aa))^2)/sum((Nmat[it,]-min(Nmat))^2)
#
K<-3
aa<-nmf(Nmat[it,]-min(Nmat),K,"brunet")
coefmap(aa)
summary(aa)
1-sum(((Nmat[it,]-min(Nmat))-fitted(aa))^2)/sum((Nmat[it,]-min(Nmat))^2)
# sparse NMF
K<-3
aa<-nmf(Nmat[it,]-min(Nmat),K,"snmf/r")
coefmap(aa)
summary(aa)
1-sum(((Nmat[it,]-min(Nmat))-fitted(aa))^2)/sum((Nmat[it,]-min(Nmat))^2)
#
K<-3
aa<-nmf(Nmat[it,]-min(Nmat),K,"snmf/l")
coefmap(aa)
summary(aa)
1-sum(((Nmat[it,]-min(Nmat))-fitted(aa))^2)/sum((Nmat[it,]-min(Nmat))^2)
########
####### can we find the digits based on the low-dim space?
kk<-kmeans(fitted(aa),10)
table(kk$cluster,Numbers[it,1])
# 
########
####### which method works best?
####
########
K<-10
VarExpl<-matrix(0,K,2)
ss<-svd(Nmat[it,]-min(Nmat))
VarExpl[,1]<-cumsum(ss$d[1:K]^2)/sum(ss$d^2)
for (kk in (1:K)) {
  aa<-nmf(Nmat[it,]-min(Nmat),kk)
  VarExpl[kk,2]<-1-sum(((Nmat[it,]-min(Nmat))-fitted(aa))^2)/sum((Nmat[it,]-min(Nmat))^2)
}
plot(seq(1,K),VarExpl[,1],ymin=c(min(VarExpl),max(VarExpl)),type="l")
lines(seq(1,K),VarExpl[,2],col=2)
####
#### SVD vs NMF
#### % explained vs interpretability
####
######
library(plot3D)
library(rgl)
library(svdvis)
ssp<-svd(puppy-mean(puppy))
svd.scree(ssp, subr=15,
          axis.title.x="Full scree plot", axis.title.y="% Var Explained")
plot3d(ssp$v[,1:3])
#### Low-rank approximations
image(puppy,col=gray.colors(256,0,1))
k<-2
image(ssp$u[,1:k]%*%diag(ssp$d[1:k])%*%t(ssp$v[,1:k]),col=gray.colors(256,0,1))
k<-5
image(ssp$u[,1:k]%*%diag(ssp$d[1:k])%*%t(ssp$v[,1:k]),col=gray.colors(256,0,1))
k<-15
image(ssp$u[,1:k]%*%diag(ssp$d[1:k])%*%t(ssp$v[,1:k]),col=gray.colors(256,0,1))
k<-25
image(ssp$u[,1:k]%*%diag(ssp$d[1:k])%*%t(ssp$v[,1:k]),col=gray.colors(256,0,1))
#### Missing values (corrupted image)
aa<-sample(seq(1,322),100,replace=T)
bb<-sample(seq(1,300),100,replace=T)
puppy2<-puppy
puppy2[aa,bb]<-NA
image(puppy2,col=gray.colors(256,0,1))
# library(softImpute)
#### Using low-rank approximations to fill in missing values
pcomp<-softImpute(puppy2,rank=10,lambda=.85)
image(pcomp$u%*%diag(pcomp$d)%*%t(pcomp$v),col=gray.colors(256,0,1))
pcomp<-softImpute(puppy2,rank=35,lambda=.9)
image(pcomp$u%*%diag(pcomp$d)%*%t(pcomp$v),col=gray.colors(256,0,1))
#### or denoising
puppy3<-puppy
puppy3<-(puppy+matrix(rnorm(322*300,sd=.12),322,300))
image(puppy3,col=gray.colors(256,0,1))
ssp<-svd(puppy3)
k<-15
image(ssp$u[,1:k]%*%diag(ssp$d[1:k])%*%t(ssp$v[,1:k]),col=gray.colors(256,0,1))
#######
### the low rank approximations of the puppy image
par(mfrow=c(4,4))
for (k in 1:16) {
  if (k>1) {
  image(ssp$u[,1:k]%*%diag(ssp$d[1:k])%*%t(ssp$v[,1:k]),col=gray.colors(256,0,1)) }
  if (k==1) {
    image(ssp$d[1]*matrix(ssp$u[,1:k],322,1)%*%matrix(ssp$v[,1:k],1,300),col=gray.colors(256,0,1)) }
}
par(mfrow=c(4,4))
for (k in 1:16) {
    image(ssp$d[k]*matrix(ssp$u[,k],322,1)%*%matrix(ssp$v[,k],1,300),col=gray.colors(256,0,1)) }
#####################
### And now NMF
aa<-nmf(puppy,16)
image(fitted(aa),col=gray.colors(256,0,1))
# coef(aa) 9*300
# basis(aa) 322*9
par(mfrow=c(4,4))
for (k in (1:16)) {
  image(matrix((basis(aa))[,k],322,1)%*%matrix((coef(aa)[k,]),1,300),col=gray.colors(256,0,1))
}
#
par(mfrow=c(4,4))
for (k in (1:16)) {
  if (k>1) {
  image((basis(aa))[,1:k]%*%(coef(aa)[1:k,]),col=gray.colors(256,0,1))
  }
  if (k==1) {
    image(matrix((basis(aa))[,k],322,1)%*%matrix((coef(aa)[k,]),1,300),col=gray.colors(256,0,1))
  }
}
###########################
#### on a big data set of digits
# this takes a while to run!
N<-5500
it<-sample(seq(1,7291),N)
K<-16
aaNbr<-nmf(Nmat[it,]-min(Nmat),K)
#
coefmap(aaNbr)
basismap(aaNbr)
###############################3
par(mfrow=c(4,4))
for (k in (1:16)) {
  image(matrix(coef(aaNbr)[k,],16,16,byrow=T),col=gray.colors(256,0,1))
}
# how good is the approximation? Look at observation u (try different ones)
# here I try the 3rd observation
u<-3
s1<-(basis(aaNbr))[u,]%*%(coef(aaNbr))
image(matrix(s1,16,16,byrow=T),col=gray.colors(256,0,1))
image(matrix((Nmat[it,])[u,],16,16,byrow=T),col=gray.colors(256,0,1))
######
svdN<-svd(Nmat[it,])
par(mfrow=c(4,4))
for (k in 1:16) {
  image(matrix(svdN$v[,k],16,16,byrow=T),col=gray.colors(256,0,1))
}
s1<-svdN$u[u,1:16]%*%diag(svdN$d[1:16])%*%t(svdN$v[,1:16])
image(matrix(s1,16,16,byrow=T),col=gray.colors(256,0,1))
image(matrix((Nmat[it,])[u,],16,16,byrow=T),col=gray.colors(256,0,1))
#
###########################
# sOM vs MDS
library(kohonen)
# 
data(crabs)
crabs2<-crabs
head(crabs)
crabs2[,1]<-as.numeric(crabs2[,1])
crabs2[,2]<-as.numeric(crabs2[,2])
pp<-prcomp(crabs2[,-3],scale=T)
CC<-as.matrix(crabs2[,-3])%*%pp$rot
plot(CC[,1],CC[,2],col=crabs2$sp,pch=crabs2$sex)
plot(CC[,3],CC[,4],col=crabs2$sp,pch=crabs2$sex)
## PCA of craps data
## It really performs poorly since the binary (categorical) data 
## is not well summarized in the SVD
ss<-som(as.matrix(crabs2[,-3]),grid=somgrid(5,5))
plot(ss,type="mapping",col=crabs2$sp,pch=crabs2$sex)
plot(ss,type="codes")
### SOMS
som_model <- som(as.matrix(crabs2[,-3]), 
                 grid=somgrid(5,5), 
                 rlen=100, 
                 alpha=c(0.05,0.01), 
                 keep.data = TRUE,
                 n.hood="circular" )
plot(som_model,type="count")
plot(som_model,type="dist.neighbours")
plot(som_model,type="codes")
#
#visualize one variable across the map
plot(som_model, type = "property", 
     property = som_model$codes[,4], 
     main=names(som_model$data)[4])
plot(som_model, type = "property", 
     property = som_model$codes[,7], 
     main=names(som_model$data)[7])
plot(som_model, type = "property", 
     property = som_model$codes[,1], 
     main=names(som_model$data)[1])
#
som_cluster <- cutree(hclust(dist(som_model$codes)), 6)
plot(som_model, type="mapping", bgcol = som_cluster, main = "Clusters") 
add.cluster.boundaries(som_model, som_cluster)
#### MDS
####
library(cluster)
dd<-daisy(crabs[,-3])
# daisy understands that some variables are categorical!
ff<-cmdscale(dd,2)
plot(ff,col=crabs2[,1],pch=crabs2[,2])
# 2d rep
#
library(plot3D)
library(rgl)
ff<-cmdscale(dd,3)
plot3d(ff,col=crabs2[,1],pch=crabs2[,2])
# In 3D we get a nice separation of species and gender!
######
## The number data (takes a while to run)
dd<-daisy(Nmat[it,])
ffNbr3<-cmdscale(dd,3)
plot3d(ffNbr3,col=as.numeric(Numbers[it,1]))
#
ffNbr2<-cmdscale(dd,2)
plot(ffNbr2,pch=as.character(Numbers[it,1]))
#
ssNbr<-som(Nmat[it,],grid=somgrid(5,5))
plot(ssNbr,type="mapping",pch=as.character(Numbers[it,1]))