#dr Jakub Paderewski
#Department of Experimental Design and Bioinformatics
#Warsaw University of Life Sciences, SGGW

AMMI.table<-function(X, err.MS=NA, err.df=NA, test="FR", rounding=c(2,2,2,3))
{
   X<-as.matrix(X)
   x.mean<-mean(X)
   X.Ie<-scale(X-x.mean,center = TRUE, scale = FALSE)
   X.Ie<-t(scale(t(X.Ie),center = TRUE, scale = FALSE))
   SVD <- La.svd(X.Ie)

   IPC.nb<-1:min(ncol(X)-1,nrow(X)-1)
   if (test=="F") {
     IPC.SS<-SVD$d[1:max(IPC.nb)]^2
     IPC.df<-nrow(X)+ncol(X)-1-2*IPC.nb
   }
   if (test=="FR") {
     IPC.SS<-rep(0,max(IPC.nb))
     IPC.df<-nrow(X)+ncol(X)-1-2*IPC.nb
     for (i in 1:max(IPC.nb))
     {
       IPC.SS[i]<-sum(SVD$d[i:max(IPC.nb)]^2)
       IPC.df[i]<-sum(IPC.df[i:max(IPC.nb)])
     }
   }
   IPC.MS<-IPC.SS/IPC.df
   IPC.F<-IPC.MS/err.MS
   IPC.pv<-pf(err.MS/IPC.MS,err.df, IPC.df)

   IPC.SS<-round(IPC.SS,rounding[1])
   IPC.MS<-round(IPC.MS,rounding[2])
   IPC.F<-round(IPC.F,rounding[3])
   IPC.pv<-round(IPC.pv,rounding[4])

   table.AMMI<-data.frame(SS=IPC.SS, df=IPC.df, MS=IPC.MS, F=IPC.F, P.v=IPC.pv)
   rownames(table.AMMI)<-paste("IPC",1:max(IPC.nb))
   return(table.AMMI)
}

mplot.AMMI<-function(X, nbPC=1, ordering="mean", cex.axis=1, File=FALSE, plot.height=7, plot.width=10, plot.res=300, plot.units="in", ...)
{
   X<-as.matrix(X)
   if (nbPC=="Full") {nbPC=min(ncol(X),nrow(X))-1}
   x.mean<-mean(X)
   X.Ie<-scale(X-x.mean,center = TRUE, scale = FALSE)
   x.col.center<-attr(X.Ie,"scaled:center")
   X.Ie<-t(scale(t(X.Ie),center = TRUE, scale = FALSE))
   x.row.center<-attr(X.Ie,"scaled:center")

   SVD <- La.svd(X.Ie)
   Gen.scores<-SVD$u[,1:nbPC]
   Env.scores<-SVD$vt[1:nbPC,]
   lambda<-diag(SVD$d[1:nbPC],nrow=nbPC)
   interaction.adj<-Gen.scores%*%lambda%*%Env.scores

   row.eff<-matrix(x.row.center,nrow(X),ncol(X))
   col.eff<-t(matrix(x.col.center,ncol(X),nrow(X)))
   Gen.yield<-interaction.adj+row.eff+col.eff+x.mean
   colnames(Gen.yield)<-colnames(X)
   rownames(Gen.yield)<-rownames(X)

   if (length(ordering)==ncol(X)) { ord.v<-ordering
   } else {
     if (ordering=="mean") { ord.v<-x.col.center }
     if (ordering=="PC1") { ord.v<-SVD$vt[1,] }
     if (ordering=="angle") { ord.v<-asin(SVD$v[2,]/ ((SVD$v[1,]^2+SVD$v[2,]^2)^0.5)) }
     if (ordering==FALSE) { ord.v<-1:ncol(X) }
   }
   ord<-order(ord.v)
 
   if (File!=FALSE)
   {
     library(grDevices)
     png(filename=File, width=plot.width, height=plot.height, res=plot.res, units=plot.units)
   }

   Environments<-1:ncol(X)
   Yield<-t(Gen.yield[,ord])
   matplot(x=Environments, y=Yield, type="b", pch=letters[1:nrow(Gen.yield)], xaxt="n", yaxt="n", col=1, lty=1, ...)
   axis(1, at=1:ncol(X), labels=LETTERS[1:ncol(Gen.yield)], lwd=1, cex.axis=cex.axis) 
   axis(2, cex.axis=cex.axis)

   if (File!=FALSE) { dev.off() }

Env.symbols=data.frame(Env=colnames(Gen.yield)[ord], ordered.by=ord.v[ord])
rownames(Env.symbols)<-LETTERS[1:ncol(Gen.yield)]
Gen.symbols=data.frame(Gen=rownames(Gen.yield))
rownames(Gen.symbols)<- letters[1:nrow(Gen.yield)]
return(list( AMMI.adjusted.means=Gen.yield, Env.symbols=Env.symbols, Gen.symbols=Gen.symbols  ))
}

nominal.yield.plot<-function(X, mean.y=TRUE, mean.col="gray", mean.lty=3, mean.lwd=3, File=FALSE, plot.height=7, plot.width=10, plot.res=300, plot.units="in", ...)
{
   X<-as.matrix(X)
   x.mean<-mean(X)
   X.Ie<-scale(X-x.mean,center = TRUE, scale = FALSE)
   x.col.center<-attr(X.Ie,"scaled:center")
   X.Ie<-t(scale(t(X.Ie),center = TRUE, scale = FALSE))
   x.row.center<-attr(X.Ie,"scaled:center")

   SVD <- La.svd(X.Ie)
   I.SS<-sum(SVD$d^2)
   lambda<-SVD$d[1]
   Gen.scores<-SVD$u[,1]
   Env.scores<-SVD$v[1,]
   interaction.adj<-lambda*(Gen.scores%*%t(Env.scores))
   ord<-order(Env.scores)

   row.eff<-matrix(x.row.center,nrow(X),ncol(X))
   col.eff<-t(matrix(x.col.center,ncol(X),nrow(X)))
   Gen.yield<-row.eff+interaction.adj
   if (mean.y!=FALSE) {Gen.yield<-Gen.yield+x.mean}
   colnames(Gen.yield)<-colnames(X)
   rownames(Gen.yield)<-rownames(X)

   if (File!=FALSE)
   {
     library(grDevices)
     png(filename=File, width=plot.width, height=plot.height, res=plot.res, units=plot.units)
   }

   Environment.range<-Env.scores[ord[c(1,NROW(ord))]]
   nominal.yield<-t(Gen.yield[,ord[c(1,NROW(ord))]])
   matplot(Environment.range, nominal.yield, type="b", pch=letters[1:nrow(Gen.yield)], lwd=2, lty=1, xaxt="n", col=1, ...)
   axis(1,at=Env.scores[ord],labels=LETTERS[1:ncol(Gen.yield)],...)
   axis(1,at=0,labels="0",...)

   if (mean.y==TRUE) {mean.y<-x.mean}
   if (!mean.y==FALSE) {abline(a=mean.y, b=0, col=mean.col, lty=mean.lty, lwd=mean.lwd)}

   if (File!=FALSE) { dev.off() }

   return(list( Env.symbols=data.frame(pch=LETTERS[1:ncol(Gen.yield)], Env=colnames(Gen.yield)[ord]), Gen.symbols=data.frame(pch=letters[1:nrow(Gen.yield)], Gen=rownames(Gen.yield)), Env.scores=data.frame(Env=colnames(Gen.yield)[ord], score=Env.scores[ord]), Gen.scores=data.frame(Gen=rownames(Gen.yield), score=Gen.scores*lambda), contained.part.of.SS=(lambda^2+sum(row.eff^2))/ (I.SS+sum(row.eff^2)) ))
}