#################################################################
# FILE: pp.ssc                                                  #
#                                                               #
# DESCRIPTION: Projection pursuit for CCA using maximization    #
#              routine                                          #
#                                                               #
# VARIABLES:  p<-dimension of the first group                   #
#             q<-dimension of the second group                  #
#             data<-matrix of the data                          #
#				inds<-indicator of the type of estimator used to   #
#                  obtain the inicial estimation of the cova-   #
#                  riance matrix                                #
#				indp<-indicator of the type of estimator of the    #
#                  projection index,                            #
#                                                               #
# OBSERVATIONS: inds,indp==1<-Classical estimator               #
#            	   inds,indp==2<-MCD estimator                     #
#           	   inds,indp==3<-M-estimator                       #
#                    indp==4<-Spearman estimator                #
#                                                               #
# LAST MODIFICATION: 3 APRIL 2002                              #
#################################################################
"pp"<-function(p,q,data,inds=1,indp=1,perc=0.75)
{
	data<-as.matrix(data)
   n<-nrow(data)

   alpha<-matrix(NA,p,p)
   beta<-matrix(NA,q,p)
   alphas<-matrix(NA,p,p)
   betas<-matrix(NA,q,p)
   ddata<-matrix(NA,n,p+q)
   tdata<-matrix(NA,n,p+q)
   sdata<-matrix(NA,n,p+q)
	base1<-matrix(NA,p,p)
	base2<-matrix(NA,q,q)
	saux<-matrix(NA,(p+q),(p+q))
	s21<-matrix(NA,q,p)
	xxtheta<-matrix(NA,p,p)
	yytheta<-matrix(NA,q,p)
	rho<-NULL
	aux1<-NULL
	aux2<-NULL
	a<-NULL
	b<-NULL
	rho<-NULL
	theta<-NULL
	xtheta<-NULL
	ytheta<-NULL
	xth<-matrix(NA,3,3)
	yth<-matrix(NA,3,3)


############ 
   dg<-diagonalise(p,q,data,inds)
	ddata<-dg$dat
	tdata<-ddata

	xxtheta<-dg$ateo
	yytheta<-dg$bteo

	for(nit in 1:(p-1))
	{
		m<-p+q-2*nit
		pp<-p-nit
		qq<-q-nit
	

		
		
######### Inicial values for alpha and beta #############
#
#		begin<-inic(pp+1,qq+1,ddata,indp)
#		xtheta<-as.vector(begin$a)
#		ytheta<-as.vector(begin$b)
#
##########################################################
##########    Alternativa ao ponto inicial ###############
#
		
		xtheta<-xxtheta[,nit]
		ytheta<-yytheta[,nit]

##########################################################


	  	theta[1:pp]<-invpolar(pp,xtheta[nit:p])
		theta[(pp+1):m]<-invpolar(qq,ytheta[nit:q])

		
#		max<-nlminb(start=theta[1:m],objective=projind,tdata=tdata[,1:(m+2)],pp=pp,qq=qq,ind=indp)
		max<-nlminb(start=theta[1:m],objective=projind,control=nlminb.control(iter.max=10000,x.tol=0.0001),tdata=tdata[,1:(m+2)],pp=pp,qq=qq,ind=indp)


		theta[1:m]<-max$parameters[1:m]

	
		aux1<-polar((p-nit+1),theta[1:(p-nit)])
		aux2<-polar((q-nit+1),theta[(p-nit+1):m])


			
		if(nit==1)
		{
	 		alphas[,1]<-(aux1)
			betas[,1]<-(aux2)
		}
		else
		{
	 		alphas[,nit]<-base1[,nit:p]%*%(aux1[1:(p-nit+1)])
			betas[,nit]<-base2[,nit:q]%*%(aux2[1:(q-nit+1)])
		}
		
		base1[,nit]<-alphas[,nit]
		base2[,nit]<-betas[,nit]	
		base1<-buildbase(p,nit,base1)
		base2<-buildbase(q,nit,base2)


		tdata[,1:(p-nit)]<-ddata[,1:p]%*%base1[1:p,(nit+1):p]
		tdata[,(p-nit+1):(p+q-2*nit)]<-ddata[,(p+1):(p+q)]%*%base2[,(nit+1):q]


		if(nit<(p-1))
		{
			xxtheta[,(nit+1)]<-t(base1)%*%dg$ateo[,(nit+1)]
			yytheta[,(nit+1)]<-t(base2)%*%dg$bteo[,(nit+1)]
		}	


	}



	base1<-buildbase(p,(p-1),base1)
	alphas[,p]<-base1[,p]
			
	if(p==q)
	{
		base2<-buildbase(q,(q-1),base2)	
		betas[,q]<-base2[,q]
	}
	else
	{
		saux<-cov.mcd(ddata,print=F,quan=n*perc)$cov
		#saux<-var(ddata)
		s21<-saux[(p+1):(p+q),1:p]
		betas[,p]<-s21%*%alphas[,p]
		betas[,p]<-betas[,p]/vecnorm(betas[,p])
	}
	

	alpha<-t(dg$c1)%*%alphas
	beta<-t(dg$c2)%*%betas



############## Calculation of the canonical correlations ############
	sdata[,1:p]<-data[,1:p]%*%alpha
	sdata[,(p+1):(p+q)]<-data[,(p+1):(p+q)]%*%beta

	for(i in 1:p)
	{
		if(indp==1)
			rho[i]<-abs( cor(cbind(sdata[,i],sdata[,(i+p)]))[1,2] )
		else
		{		
		#	s<- mesthub(cbind(sdata[,i],sdata[,(i+p)]),delta=0.05)$cov	
		#	rho[i]<-abs(s[1,2]/sqrt(s[1,1]*s[2,2]))		
		#	s<- cor(cbind(sdata[,i],sdata[,(i+p)]))
		perc<-0.75
		s<-cov.mcd(cbind(sdata[,i],sdata[,(i+p)]),print=F,cor=T,quan=floor(n*perc))
		rho[i]<-s$cor[1,2]
			
		}
	}


#		rho[i]<-abs(cov.mcd(cbind(sdata[,i],sdata[,(i+p)]),print=F,cor=T,quan=0.75*n)$cor[1,2])



	list(a=alpha,b=beta,rho=rho)

}





#################################################################
# SUBROUTINE: inic                                              #
#                                                               #
# DESCRIPTION: calculates the inicial value for alpha and beta  #
#                                                               #
# VARIABLES:  p<-number of variables of the first group of      #
#                variables                                      #
#             q<-number of variables of the second group of     #
#                variables                                      #
#             data<- data                                       #
#				 ind<-type of estimator used for the projection    #
#                  index                                        #
#################################################################
"inic"<-function(p,q,data,ind=1)
{
	n<-nrow(data)
	
	dat.x<-matrix(NA,n,p)
	dat.y<-matrix(NA,n,q)
	cdat.x<-matrix(NA,n,p)
	cdat.y<-matrix(NA,n,q)
	mx<-matrix(NA,n,n)
	vpos<-NULL
	loc<-NULL
	maxi<-NULL
	aux1<-NULL
	aux2<-NULL
	
	
	nind<-2
	
	
	dat.x <- data[1:n, 1:p]
	dat.y <- data[1:n, (p + 1):(p + q)]
		
	#loc<-cov.mcd(data,quan=n*perc)$center
	loc<-apply(data,2,median)
	cdat.x<-scale(dat.x, center=loc[1:p], scale=F)
	cdat.y<-scale(dat.y, center=loc[(1+p):(p+q)], scale=F)
	

	for(i in 1:n)
	{
		cdat.x[i,]<-cdat.x[i,]/vecnorm(cdat.x[i,])
		cdat.y[i,]<-cdat.y[i,]/vecnorm(cdat.y[i,])
	}
	
		
	
	if(nind==1)
	{
		for(i in 1:n)
			for(j in 1:n)
			{
				aux1<-dat.x%*%cdat.x[i,]
				aux2<-dat.y%*%cdat.y[j,]
				mx[i,j]<-IP(aux1,aux2,ind)
			}
		vpos<-maxindc(mx)
	
		alpha<-as.vector(cdat.x[vpos[1],])
	 	beta<-as.vector(cdat.y[vpos[2],])
	}
	else
	{	
		pos<-1
		aux1<-dat.x%*%cdat.x[1,]
		aux2<-dat.y%*%cdat.y[1,]
		mx<-IP(aux1,aux2,ind)
		
		for(i in 2:n) 
		{
			aux1<-dat.x%*%cdat.x[i,]
			aux2<-dat.y%*%cdat.y[i,]
			maxi[i]<-IP(aux1,aux2,ind)
			if(maxi[i]>mx)
			{
				pos<-i
				mx<-maxi[i]
			}
		}
   		alpha<-as.vector(cdat.x[pos,])
   		beta<-as.vector(cdat.y[pos,])
	}
	
	list(a = alpha , b = beta)
}




#################################################################
# SUBROUTINE: IP                                                #
#                                                               #
# DESCRIPTION: Calculates the projection index                  #
#                                                               #
# VARIABLES:  x<-first set of data                              #
#             y<-second set of data                             #
#				 ind<-type of estimator used for the projection    #
#                  index                                        #
#################################################################
"IP"<-function(x,y,ind=1,perc=0.75)
{
	s<-matrix(NA,2,2)
	n<-nrow(x)
	xy<-matrix(NA,n,2)
	
	xy<-cbind(x,y)

	if(ind==1)
		proj<-var(xy)[1,2]
	else if(ind==2)
	{
		ss<-cov.mcd(xy,print=F,cor=T,quan=floor(n*perc))
		proj<-ss$cor[1,2]
	}
	else if(ind==3)
	{
		s<-mesthub(xy,delta=0.05)$cov
		proj<-s[1,2]/sqrt(s[1,1]*s[2,2])
	}
	else
#		proj<-cor(apply(xy,2,rank))[1,2]
		proj<-2*sin(pi*cor(apply(xy,2,rank))[1,2]/6)
		
	return(proj)
}





##################################################################
# SUBROUTINE: diagonalise                                        #
#                                                                #
# DESCRIPTION: Given the data their covariance matrix is estima- #
#              ted according with ind.                           #
#				  The covariance matrix of each group are facto-    #
#              rized in the form: Sigma_ii=SS' in order to obtain#
#              an equivalent problem expressed in euclidian      #
#              spaces                                            #
#             															  #
# VARIABLES:  p<-number of variables of the first group          #
#             q<-number of variables of the second group         #
#             data<-original data                                #
#             ind<-index of the type of estimator to be used     # 
##################################################################
"diagonalise"<-function(p,q,data,ind=1,perc=0.75)
{
	n<-nrow(data)
	ddata<-matrix(NA,n,p+q)
	sigma<-matrix(NA,p+q,p+q)
	ateo<-matrix(NA,p,p)
	bteo<-matrix(NA,q,q)
	
	
######## Inicial estimate of the covariance matrix ##########
	if(ind==1) 
		sigma<-var(data) 
	else if(ind==2) 
	{		
		ss<-cov.mcd(data,print=F,cor=T,quan=n*perc)
		sigma<-ss$cov
	}
	else
		sigma<-mesthub(data,0.05)$cov



########### Calculation of c1 and ic1 #############
	eigsysx <- eigen(sigma[1:p,1:p])
	vectx <- eigsysx$vectors
	valx<-eigsysx$values
	ia1<-diag(sqrt(valx))
	a1<-diag(1/sqrt(valx))
	c1<-a1%*%t(vectx)
	ic1<-ia1%*%t(vectx)
	
	ddata[,1:p]<-data[,1:p]%*%t(c1)
	
	xtheta<-diag(p)
	ateo<-ic1%*%xtheta

	
	for(i in 1:p)
		ateo[,i]<-ateo[,i]/vecnorm(ateo[,i],p=2)
		

	
########### Calculation of c2 and ic2 #############	
	eigsysy <- eigen(sigma[(1+p):(p+q),(1+p):(p+q)])
	vecty <- eigsysy$vectors
	valy<-eigsysy$values
	ia2<-diag(sqrt(valy))
	a2<-diag(1/sqrt(valy))
	c2<-a2%*%t(vecty)
	ic2<-ia2%*%t(vecty)
	

	ddata[,(1+p):(p+q)]<-data[,(1+p):(p+q)]%*%t(c2)
	
	ytheta<-diag(q)
	bteo<-ic2%*%ytheta
	for(i in 1:q)
		bteo[,i]<-bteo[,i]/vecnorm(bteo[,i],p=2)



		
   list(dat=ddata,c1=c1,ic1=ic1,c2=c2,ic2=ic2,ateo=ateo,bteo=bteo)
}



#################################################################
# SUBROUTINE: invpolar                                          #
#                                                               #
# DESCRIPTION: Given a vector x, transformes it to polar coorde-#
#              nates                                            #
#                                                               #
# VARIABLES:  k<-dimension of the vector x                      #
#             x<-vector to be transform                         #
#################################################################
"invpolar"<-function(k,x)
{
	theta<-NULL
	
	if(k==1)
		theta[1]<-atan(x[2]/x[1])
	else if(k==2)
	{
		theta[1]<-atan(x[2]/x[1])
		theta[2]<-atan(x[2]/(x[3]*sin(theta[1])) )
	}
	else
	{
		theta[1]<-atan(x[2]/x[1])
		theta[2]<-atan(x[2]/(x[3]*sin(theta[1])) )
		for(i in 3:k)
		{
			theta[i]<-atan(x[i]/(x[i+1]*cos(theta[i-1])))
		}
	
	}
	return(theta)
}



#################################################################
# SUBROUTINE: polar                                             #
#                                                               #
# DESCRIPTION: Given a vector in polar coordenates transformes  #
#              it to the original space                         #
#                                                               #
# VARIABLES:  k<-dimension of theta+1                           #
#             theta<-vector in polar coordenates                #
#################################################################
"polar"<-function(k,theta)
{
	x<-NULL
	aux<-NULL
	
	if(k==2)
	{
		x[1]<-cos(theta[1])
		x[2]<-sin(theta[1])	
	}
	else 
	{
		
		aux<-sin(theta[k-1])*polar(k-1,theta[1:(k-2)])
		x[1:(k-1)]<-aux[1:(k-1)]
		x[k]<-cos(theta[(k-1)])
	}

	return(x[1:k])
}


#################################################################
# SUBROUTINE: projind                                           #
#                                                               #
# DESCRIPTION: Objective function to be maximized by nlminb, in #
#              order to maximize the projection index according #
#              with the index: ind                              #
#                                                               #
# VARIABLES:  theta<-vector where the projection index should be#
#                    evaluated, written in polar coordenates    # 
#             tdata<-data where the projection should be eva-   #
#                    luated                                     #
#             pp<-dimension in use for the first group of       #
#                 variables                                     #
#             qq<-dimension in use for the second group of      #
#                 variables 
#             ind<-indicator of the estimator for the pro-      #
#                  jection index
##################################################################
"projind"<-function(theta,tdata,pp,qq,ind,perc=0.75)
{
	n<-nrow(tdata)
	z<-matrix(NA,n,2)
	s<-matrix(NA,2,2)
	alpha<-NULL
	beta<-NULL
	theta<-as.vector(theta)
	
	n<-nrow(tdata)
	m<-pp+qq
	alpha<-polar((pp+1),theta[1:pp])
	beta<-polar((qq+1),theta[(pp+1):m])
	

	z[,1]<-tdata[,1:(pp+1)]%*%alpha
	z[,2]<-tdata[,(pp+2):(m+2)]%*%beta

	
	if(ind==1)
	{
		s<-cor(z)
		fc<--s[1,2]
	}
	else if (ind==2)
	{
		s<-cov.mcd(z,print=F,cor=T,quan=n*perc)
		fc<--s$cor[1,2]
	}
	else if(ind==3)
	{
		s<-mesthub(z,0.05)$cov
		fc<--s[1,2]/sqrt(s[1,1]*s[2,2])		
	}
	else
		fc<--cor(apply(z,2,rank))[1,2]
#		fc<--2*sin(pi*cor(apply(z,2,rank))[1,2]/6)

	return(fc)
	
	
}


#################################################################
# SUBROUTINE: buildbase                                         #
#                                                               #
# DESCRIPTION: Given a set of k independent vector builds a base#
#              in a space of dimension p having those as the    #
#              first k elements                                 #
#                                                               #
# VARIABLES:  p<-dimension of the space                         #
#             k<-number of independent vectors                  #
#             base<-independent vectors                         #
#################################################################
"buildbase"<-function(p,k,base)
{
	base<-as.matrix(base)
	aux<-matrix(NA,p,p)
	id<-matrix(NA,p,p)
	
	id<-diag(p)
	aux[1:p,1:p]<-id[1:p,1:p]-base[1:p,1:k]%*%t(base[1:p,1:k])
	
	naux<-0
	i<-1
	
	while((naux!=(p-k))&&(i<=p))
	{
		x<-aux%*%id[,i]
		if(all(abs(x)<0.00001)) 
			i<-i+1
		else
		{
			naux<-naux+1
			base[,(k+naux)]<-id[,i]
			base<-GSchmidt((k+naux),base)
			i<-i+1
			aux[1:p,1:p]<-id[1:p,1:p]-base[1:p,1:(k+naux)]%*%t(base[1:p,1:(k+naux)])
		}
	}
	return(base)
}



#################################################################
# SUBROUTINE: GSchmidt                                          #
#                                                               #
# DESCRIPTION: Performs orthogopmalization of Gram Schmidt,     #
#              return orthornoamal vectors.                     #
#                                                               #
# VARIABLES:  k<-number of independent vectors                  #
#             base<-independent vectors                         #
#################################################################
"GSchmidt"<-function(k,base)
{
	p<-nrow(base)
	new<-matrix(NA,p,p)
	
	new<-base
	
	for(i in 2:k)
	{
		for(j in 1:(i-1))
			new[,i]<-new[,i]-(t(base[,i])%*%new[,j]/vecnorm(new[,j])^2)*base[,j]
		new[,i]<-new[,i]/vecnorm(new[,i])
	}
	
	return(new)	
	
}