## Splus functions for depicting ternary diagrams
## Function "two.dim" is used within "ternary.graph"
## Function "ternary.graph" depicts compositional data in a ternary diagram in a pdf file.

## Example 1
   x1 <- seq(0,1,.005)
   b0 <- 2
   b1 <- 1
   linear.regression <- exp(b0 + b1* log(x1/(1-x1)))
   x2 <- (linear.regression/(1+linear.regression))*(1-x1)
   x3 <- 1 - x1 - x2
   X.on.regression.line <- cbind(x1,x2,x3)
   ternary.graph(X.on.regression.line,filename="example1.pdf",Xlinetype="l")

## Example 2
   X <- matrix(c(.2,.2,.6,.3,.3,.4,.1,.5,.4),3,3,byrow=T)
   ternary.graph(X,filename="example2.pdf", Xpch=15, budget.label=c("Top","Left","Right"))

## Example 3
   ternary.graph(X,Y=X.on.regression.line,filename="example3.pdf", Xpch=15, budget.label=c("Top","Left","Right"))

## If there are more than three components, components 4, ..., D are concatenated.
## If there are less than three components, the function does nothing.
## X and filename are compulsory
## "Y" is a second compositional data matrix, for example, to draw a regression line in a data plot (Example 3).
## "grid" is the scale of the grid.
## "Xlinetype" is the type for data X
## "Ylinetype" is the type for data Y
## Xpch is the pch for data X
## budget.labels are the budget labels

"two.dim"<- function(x) {
    cbind(x[,2]*sqrt(4/3) + sqrt(x[,1]^2 * (1/3)), x[,1])
}

"ternary.graph"<- function(X, # data
         filename,
         Y=matrix(c(1,0,0),1,3), # estimated data
         grid = seq(0,1,.1),
         Xlinetype="p",
         Ylinetype="l",
         Xpch=1,
         budget.label=c("1","2","3"))
{
    ncol.X <- ncol(X)
    ncol.Y <- ncol(Y)
    if(ncol.X > 2){
        n.grid <- length(grid)
        # Join columns 3:ncol.X
        if(ncol.X > 3){
            X.3 <- cbind(X[,1:2],X[,3:ncol.X] %*% rep(1,ncol.X-2))
        }else {
          X.3 <- X
        }
        # Join columns 3:ncol.Y
        if(ncol.Y > 3){
            Y.3 <- cbind(Y[,1:2],Y[,3:ncol.Y] %*% rep(1,ncol.Y-2))
        }else{
          Y.3 <- Y
        }

        # Construct proper budgets
        X.3 <- X.3 / rep(X.3 %*% rep(1,3),3)
        Y.3 <- Y.3 / rep(Y.3 %*% rep(1,3),3)

       budgets_two.dim(diag(3))

       data <- two.dim(X.3)
       parameters <- two.dim(Y.3)
       scale <- 0.5
       pdf.graph(filename,horizontal=T)
       par(mai=c(.5,.5,.5,.5),pin=c(7,7))
     # plot budgets
       plot(data[,1],data[,2],
        xlim=c(0,sqrt(4/3)),ylim=c(0,1),
        xlab="",
        ylab="",
        axes=F,
        type=Xlinetype,
        pch=Xpch,
        lwd=3)
       if(nrow(Y) > 1)lines(parameters[,1],parameters[,2], type=Ylinetype, lwd=3)
     # plot outer triangle
       segments(budgets[1,1],budgets[1,2],budgets[2,1],budgets[2,2],lwd=2)
       segments(budgets[2,1],budgets[2,2],budgets[3,1],budgets[3,2],lwd=2)
       segments(budgets[3,1],budgets[3,2],budgets[1,1],budgets[1,2],lwd=2)
     # plot grid
        budget2.scale <- two.dim(cbind(rep(0,n.grid),1-grid,grid))
        budget1.scale <- two.dim(cbind(grid,rep(0,n.grid),1-grid))
        budget3.scale <- two.dim(cbind(,1-grid,grid,rep(0,n.grid)))
       tmp_cbind(two.dim(cbind(rep(0,n.grid),1-grid,grid)),two.dim(cbind(grid,1-grid,rep(0,n.grid))))
       segments(tmp[,1],tmp[,2],tmp[,3],tmp[,4],lty=2,lwd=1)
       tmp_cbind(two.dim(cbind(rep(0,n.grid),1-grid,grid)),two.dim(cbind(1-grid,rep(0,n.grid),grid)))
       segments(tmp[,1],tmp[,2],tmp[,3],tmp[,4],lty=2,lwd=1)
       tmp_cbind(two.dim(cbind(grid,1-grid,rep(0,n.grid))),two.dim(cbind(grid,rep(0,n.grid),1-grid)))
       segments(tmp[,1],tmp[,2],tmp[,3],tmp[,4],lty=2,lwd=1)
     # plot budget labels
       label.space <- .1
       text(0.5*sqrt(4/3),1+label.space, budget.label[1], cex=1.5)
       text(-label.space,0, budget.label[2], cex=1.5)
       text(sqrt(4/3)+label.space,0, budget.label[3], cex=1.5)
     # plot scales
       if(n.grid > 2){
         for(i in 2:(n.grid-1)){
          text(budget1.scale[i,1]-0.5*label.space,budget1.scale[i,2],as.character(grid[i]),cex=1.5)
          text(budget2.scale[i,1],budget2.scale[i,2]-0.5*label.space,as.character(grid[i]),cex=1.5)
          text(budget3.scale[i,1]+0.5*label.space,budget3.scale[i,2],as.character(grid[i]),cex=1.5)
         }
       }

      # arrows(budget1.scale[6,1]- 1.5*label.space, budget1.scale[6,2],budget1.scale[7,1]-1.5*label.space, budget1.scale[7,2],lwd=4)
      # arrows(budget2.scale[6,1],budget2.scale[6,2]-1.5*label.space,budget2.scale[7,1],budget2.scale[7,2]-1.5*label.space,lwd=4)
      # arrows(budget3.scale[6,1]+ 1.5*label.space, budget3.scale[6,2], budget3.scale[7,1]+ 1.5*label.space, budget3.scale[7,2],lwd=4)
       dev.off()
    }
}