lfda

## ruiling@ 07/08/2012
## input:
##		x: 		sample data matrix, d×n
##					d is the dimension of original sample
##					n is the number of samples
##		y: 		class label vector, n×1
##		r: 		reduced dimension
##		metric：type of metric in the embedding space
##					'weighted'			------ weighted eigenvectors
##					'orthonormalized'	------ orthonormalized
##					'plain'				------ raw eigenvectors
## 		knn: 	parameter used in local scaling method(default:7)
##
## ouput:
##		T: 		d×r transformation matrix (Z=T^(-1)*x)
##		Z: 		r×n matrix of dimensionality reduce samples
##

##------------------------------------------------------
##	test data
##------------------------------------------------------

#   x=matrix(c(1:20),4,5)
#	x=x[1:2,]
#	y=c(1,2,2,1,2)
#	r=2
#	metric='weighted'
#	k=1

#	setwd("C:/Users/ruiling/Desktop/mylfda")
	
##-------------------------------------------------------------
##  connect to Matlab Server and call Matlab function
##	step1 load R.matlab
##	step2 Start Matlab
##	step3 Create a Matlab clent object used to communicate with Matlab
##  step4 connect to the Matlab server
##---------------------------------------------------------------
#	library(R.matlab)
#	Matlab$startServer()
#	matlab <- Matlab(host="localhost",port=9999)
#	Isopen<-open(matlab)

 
##------------------------------------------------------------
## repeat a matrix like the MATLAB grammar
##------------------------------------------------------------


repmat<-function(A,N,M){
	kronecker(matrix(1,N,M),A)
}




lfda<-function(x,y,r,metric=c('weighted', 'orthonormalized', 'plain'),knn=7){

##	metric='weighted'
	n=ncol(x)
	d=nrow(x)
	if(is.null(r))	r=d
	metric=match.arg(metric)
	
	tSb=mat.or.vec(d,d)		
	tSw=mat.or.vec(d,d)
	
	for(c in unique(y)){
	
		xc=x[,y==c]			# d×nc
		nc=dim(xc)[2]		
	
	##------- Define classwise affinity matrix A --------
	
		A=mat.or.vec(nc,nc)
		xc2=t(as.matrix(colSums(xc^2)))
		distance2=repmat(xc2,nc,1)+repmat(t(xc2),1,nc)-2 * t(xc) %*% xc
		sorted=apply(distance2,2,sort)
		index=apply(distance2,2,order)
		KNNdist2=t(as.matrix(sorted[knn+1, ]))
		sigma=sqrt(KNNdist2)
		localscale=t(sigma) %*% sigma
		flag=(localscale != 0)
		A[flag]=exp(-distance2[flag]/localscale[flag])
	
	##------ Within-class & between-class scatter matrix -----
		
		xc1=as.matrix(rowSums(xc))
		G=xc %*% (repmat(as.matrix(colSums(A)),1,d)*t(xc))-xc %*% A %*% t(xc)
		tSb=tSb+G/n+xc%*%t(xc)*(1-nc/n)+xc1%*%t(xc1)/n
		tSw=tSw+G/nc
	}
	
	x1=as.matrix(rowSums(x))
	tSb=tSb-x1%*%t(x1)/n-tSw
	
	tSb=(tSb+t(tSb))/2
	tSw=(tSw+t(tSw))/2
	
	setVariable(matlab,tSb=tSb,tSw=tSw,r=r)
	
	if(r==d){
		
		evaluate(matlab,"[eigvec,eigval]=eig(tSb,tSw);")
		eigVec=getVariable(matlab,"eigvec")$eigvec
		eigVal_matrix=getVariable(matlab,"eigval")$eigval
		
	}else{
		
		evaluate(matlab,"[eigvec,eigval]=eigs(tSb,tSw,r,'la');")
		eigVec=getVariable(matlab,"eigvec")$eigvec
		eigVal_matrix=getVariable(matlab,"eigval")$eigval
	}
	
	eigval=diag(eigVal_matrix)
	
	sort_eigval=apply(t(eigval),1,sort)
	sort_eigval_index=apply(t(eigval),1,order)
	 
	
	T0 = eigVec[,sort_eigval_index[1:r]]
	T0=as.matrix(T0)
	
	T=switch(metric,
		weighted=T0*repmat(t(sqrt(eigval)),d,1),
		orthonormalized = qr.Q(qr(T0)),
		plain=T0
	)
	
	Z=t(T)%*%x
	
	
	return(list("T" = T, "Z" = Z))
}

##close(matlab)