################################################################################
# Fucntion: dlamb
# Description: function which compute generalized average
# Input:  NXM matrix with genes and conditions; degree for average;
# degree=0, -100, 100 represent limits when degree->0, -inf and +inf, respectively
# Output: MxM distance matrix
# Need: library('e1071')
################################################################################

dlamb<-function(d, degree)
{
n<-dim(d)
dm<-as.matrix(d)
dt<-as.matrix(t(d))
dmn<-dt%*%dm
da<-dmn
da[,]<-0

for(i in 1:n[2])
for(j in 1:n[2])
 {
 if(degree!=0&&degree!=100&&degree!=-100)
  tmp<-((dmn[i,i]^degree + dmn[j,j]^degree)/2)^(1/degree)
 else
 {
  if(degree==0)
    tmp<-(dmn[i,i]*dmn[j,j])^(1/2)
  if(degree==100)
    tmp<-max(dmn[i,i], dmn[j,j])
  if(degree==-100)
    tmp<-min(dmn[i,i],dmn[j,j])
  }
 da[i,j]<-dmn[i,j]/tmp
  }
return (da)
}
################################################################################
# we start here:

 degree<-readline("Please input GA distance exponent (lambda):")
 degree<-as.numeric(degree)

 finp<-readline("Please indicate the input file name for your data matrix (tab.delimited \'\\t\'):")
 d<-read.table(finp,row.names=1,header=TRUE)
foutp<-paste(finp,".lambda",sep="")
 dat<-t(d)
 n<-dim(dat)
 n1<-n[1]
 dat<-1-dlamb(d,degree)
#output lamda distance
cat("Distance matrix for your lamda = ",degree,"\n",file=foutp,append=TRUE)
write.table(formatC(dat,format="f",digit=3),file=foutp,sep='\t',quote=FALSE)

m<-2*degree+4 # 3 boundary cases and correlation
dv<-matrix(0,m,n1*(n1-1)/2)
ind<--degree:degree
for(i in 1:(m-3))
{  dat<-1-dlamb(d,ind[i])
   dv[i,]<-dat[lower.tri(dat)]
  }

library('e1071')
dst<-matrix(0,m,5)
for(i in 1:(m-3))
dst[i,]<-cbind(mean(dv[i,]),var(dv[i,]),median(dv[i,]),kurtosis(dv[i,]),skewness(dv[i,]))

dcor<-1-cor(d)
d_inf<-1-dlamb(d,-100)
dinf<-1-dlamb(d,100)
dcor<-as.vector(dcor[lower.tri(dcor)])
d_inf<-as.vector(d_inf[lower.tri(d_inf)])
dinf<-as.vector(dinf[lower.tri(dinf)])

dst[m-2,]<-rbind(mean(dcor),var(dcor),median(dcor),kurtosis(dcor),skewness(dcor))
dst[m-1,]<-rbind(mean(dinf),var(dinf),median(dinf),kurtosis(dinf),skewness(dinf))
dst[m,]<-rbind(mean(d_inf),var(d_inf),median(d_inf),kurtosis(d_inf),skewness(d_inf))
colnames(dst)<-cbind("Dist\tMean","Var","Median","Kurtosis","Skewness")
row.names(dst)<-c(-degree:degree,"dcor","dinf","d_inf")
foutp<-paste(finp,".stat",sep="")
write.table(formatC(dst,format="f",digit=3),file=foutp,sep='\t',quote=FALSE)
################################################################################
#here we select the 'best performing' distance measure: maximum kurtosis, skew

ismax<-which.max(abs(dst[,5]))#index of abs. max. skewness
ikmax<-which.max(abs(dst[,4]))#index of abs. max. kurtosis

foutpS<-paste(finp,".maxSk",sep="")
foutpK<-paste(finp,".maxKu",sep="")

#print dist with max abs. sk
if(ismax!=(m-2)&&ismax!=(m-1)&&ismax!=m)
{ #ismax-(degree+1)
  cat("Distance with abs. max. skeweness: lambda = ",ind[ismax],"\n",file=foutpS,append=TRUE)
  dat<-1-dlamb(d,ind[ismax])
  write.table(formatC(dat,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
  }
if(ismax==(m-2))
{ cat("Distance with abs. max. skeweness: d = 1-corr","\n",file=foutpS,append=TRUE)
  dcor<-1-cor(d)
  write.table(formatC(dcor,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
  }
if(ismax==m)
{ cat("Distance with abs. max. skeweness: lambda = -infinity (Simpson similarity index)","\n",file=foutpS,append=TRUE)
  d_inf<-1-dlamb(d,-100)
  write.table(formatC(d_inf,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
}
if(ismax==(m-1))
{cat("Distance with max. skeweness: lambda = +infinity ","\n",file=foutpS,append=TRUE)
 dinf<-1-dlamb(d,100)
 write.table(formatC(dinf,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
}
#and the same for kurtosis
if(ikmax!=(m-2)&&ikmax!=(m-1)&&ikmax!=m)
{ cat("Distance with abs. max. kurtosis: lambda = ", ind[ikmax],"\n",file=foutpK,append=TRUE)
  dat<-1-dlamb(d,ind[ikmax])
  write.table(formatC(dat,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
  }
if(ikmax==(m-2))
{ cat("Distance with abs. max. kurtosis: d = 1-corr","\n",file=foutpK,append=TRUE)
  dcor<-1-cor(d)
  write.table(formatC(dcor,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
  }
if(ikmax==m)
{cat("Distance with abs. max. kurtosis: lambda = -infinity (Simpson similarity index)","\n",file=foutpK,append=TRUE)
 d_inf<-1-dlamb(d,-100)
 write.table(formatC(d_inf,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
}
if(ikmax==(m-1))
{ cat("Distance with abs. max. kurtosis: lambda = +infinity ","\n",file=foutpK,append=TRUE)
  dinf<-1-dlamb(d,100)
  write.table(formatC(dinf,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
}

################################################################################
#here we select the 'best performing' distance measure: minimum kurtosis, skew

ismin<-which.min(abs(dst[,5]))#index of abs. min. skewness
ikmin<-which.min(abs(dst[,4]))#index of abs. min. kurtosis

foutpS<-paste(finp,".minSk",sep="")
foutpK<-paste(finp,".minKu",sep="")

#print dist with max abs. sk
if(ismin!=(m-2)&&ismin!=(m-1)&&ismin!=m)
{ cat("Distance with abs. min. skeweness: lambda = ", ind[ismin],"\n",file=foutpS,append=TRUE)
  dat<-1-dlamb(d,ind[ismin])
  write.table(formatC(dat,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
  }
if(ismin==(m-2))
{ cat("Distance with abs. min. skeweness: d = 1-corr","\n",file=foutpS,append=TRUE)
  dcor<-1-cor(d)
  write.table(formatC(dcor,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
  }
if(ismin==m)
{ cat("Distance with abs. min. skeweness: lambda = -infinity (Simpson similarity index)","\n",file=foutpS,append=TRUE)
  d_inf<-1-dlamb(d,-100)
  write.table(formatC(d_inf,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
}
if(ismin==(m-1))
{cat("Distance with min. skeweness: lambda = +infinity ","\n",file=foutpS,append=TRUE)
 dinf<-1-dlamb(d,100)
 write.table(formatC(dinf,format="f",digit=3),file=foutpS,quote=FALSE,sep='\t',append=TRUE)
}
#and the same for kurtosis
if(ikmin!=(m-2)&&ikmin!=(m-1)&&ikmin!=m)
{ cat("Distance with abs. min. kurtosis: lambda = ", ind[ikmin],"\n",file=foutpK,append=TRUE)
  dat<-1-dlamb(d,ind[ikmin])
  write.table(formatC(dat,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
  }
if(ikmin==(m-2))
{ cat("Distance with abs. min. kurtosis: d = 1-corr","\n",file=foutpK,append=TRUE)
  dcor<-1-cor(d)
  write.table(formatC(dcor,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
  }
if(ikmin==m)
{cat("Distance with abs. min. kurtosis: lambda = -infinity (Simpson similarity index)","\n",file=foutpK,append=TRUE)
 d_inf<-1-dlamb(d,-100)
 write.table(formatC(d_inf,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
}
if(ikmax==(m-1))
{ cat("Distance with abs. min. kurtosis: lambda = +infinity ","\n",file=foutpK,append=TRUE)
  dinf<-1-dlamb(d,100)
  write.table(formatC(dinf,format="f",digit=3),file=foutpK,quote=FALSE,sep='\t',append=TRUE)
}