#################
##  No one should alter the following, exept T.Yano, T.Omori, Y.Fujimoto.
##  If you use this script outside the Quantitative Analysis course, 
##  CIS Doshisha University, you should cite 
##  "Diagram by the script 'PCA.diagram', Quantitative Analysis course, 
##  CIS Doshisha University."
##  The following scripts are based on fa.diagram and dia.arrow.

PCA.diagram <- function(
	data,
	ncp=5,
	Phi = NULL, 
	sort = TRUE, 
	labels = NULL, 
    cut = 0.3, 
	simple = TRUE, 
	errors = FALSE, 
	digits = 2, 
	e.size = 0.05, 
    rsize = 0.2,
	side = 2, 
	main = NULL, 
	cex = 1, 
	col2="red",
    xlab="", 
	ylab="", 
	sub="", 
	graph=FALSE, 
	protate=NULL,
	pm="Facto", 
	n.obs=NULL, 
	max.cp=5,
    fa=FALSE, 
	nfactors=NULL, 
	fm="minres", 
    rotate="oblimin",
    new.plot=TRUE, 
	addP=NULL, 
	addF=NULL, 
	lmove=0, 
	...)
{
## ver. 0.97   2012.7.1
## version 0.95    2011.8.12
## version 0.96    2011.8.31     output rotated.scores

ver <- "0.97"

  require(FactoMineR)
  require(psych)

  data.name <- deparse(substitute(data))
  if (is.null(nfactors)) nfactors <- ncp

  if (is.null(main)) 
    main <- data.name

 dia.arrow2 <-
 function (from, to, labels = NULL, scale = 1, cex = 1, lmove=0,
        digits, ...) 
 {
    radius1 <- radius2 <- 0
    if (is.list(from)) {
        if (!is.null(from$radius)) {
            radius1 <- from$radius
            radius2 <- 0
            from <- from$center
        }
    }
    if (is.list(to)) {
        if (!is.null(to$radius)) {
            radius2 <- to$radius
            to <- to$center
        }
    }
    theta <- atan((from[2] - to[2])/(from[1] - to[1]))
    x <- (from[1] + to[1])/2
    y <- (from[2] + to[2])/2
    if (is.null(labels)) {
        h.size <- 0
    }
    else {
        h.size <- nchar(labels) * cex * 0.4
    }
    if (is.null(labels)) {
        v.size <- 0
    }
    else {
        v.size <-  cex * 0.2
    }
    if (from[1] < to[1]) {
        h.size <- -h.size
        radius1 <- -radius1
        radius2 <- -radius2
    }
    x0 <- from[1] - cos(theta) * radius1
    y0 <- from[2] - sin(theta) * radius1
    xe <- to[1] + cos(theta) * radius2
    ye <- to[2] + sin(theta) * radius2
    xm <- x + lmove * cos(theta)
    ym <- y + lmove * sin(theta)
    xr <- xm + h.size * cos(theta) *v.size
    xl <- xm - h.size * cos(theta) *v.size
    yr <- ym + v.size * sin(theta) *h.size
    yl <- ym - v.size * sin(theta) *h.size
    if (!is.null(labels)) {
        text(xm, ym, labels, cex = cex, ...)
    }
    arrows(x0, y0, xr, yr, length = 0, angle = 30, code = 2, 
        ...)
    arrows(xl, yl, xe, ye, length = 0.1 * scale, angle = 30, 
        code = 2, ...)

}  # end dia.arrow2


 pca.diagram <-
 function (fa.results,prop, Phi, sort, labels, 
    cut, simple, errors, digits, e.size, 
    rsize, side, main, cex, col2,protate,
    xlab, ylab, sub, data.name, max.cp, 
    fa, nfacors, fm, rotate, new.plot, n.obs, addF, addP, lmove, ...) 
{
    col <- c("black", col2)

#    if (is.null(cex)) 
#        cex <- 1
    if (sort) 
        fa.results <- fa.sort(fa.results)
    if ((!is.matrix(fa.results)) && (!is.data.frame(fa.results))) {
        factors <- as.matrix(fa.results$loadings)
        if (!is.null(fa.results$Phi)) {
            Phi <- fa.results$Phi
        }
        else {
            if (!is.null(fa.results$cor)) {
                Phi <- fa.results$cor
            }
        }
    }
    else {
        factors <- fa.results
    }
    nvar <- dim(factors)[1]  ########
    if (is.null(nvar)) {
        nvar <- length(factors)
        num.factors <- 1
    }
    else {
        num.factors <- dim(factors)[2]
    }
    nvar <- dim(factors)[1]
    e.size = e.size * 16/nvar
    if (is.null(nvar)) {
        nvar <- length(factors)
        num.factors <- 1
    }
    else {
        num.factors <- dim(factors)[2]
    }
    if (is.null(rownames(factors))) {
        rownames(factors) <- paste("V", 1:nvar, sep = "")
    }
    if (is.null(colnames(factors))) {
        colnames(factors) <- paste("F", 1:num.factors, sep = "")
    }
    var.rect <- list()
    fact.rect <- list()
    max.len <- max(nchar(rownames(factors))) * rsize
    x.max <- max((nvar + 1), 6)
    limx = c(-max.len*1.1, x.max)
    limy = c(1, nvar)
    unique <- 1 - rowSums(factors^2)

  if (!fa) {
       if (pm=="Facto" & is.null(n.obs)) main <- paste("PCA diagaram:" ,main)
       else main <- paste("PCA diagaram:" ,main,", ",protate)
       xlab <- paste(xlab,"PCA loadings, ver.",ver,". left:residuals, right:cumulative proportions.")
  } else {
       main <- paste(main, "Path diagaram" ,fm,rotate,
               ", nf =", nfactors)
       xlab <- paste(xlab,"FA path diagram. ver.",ver,". left:uniquenesses, right:cum.prop.")
  }

  if (new.plot) dev.new()
    plot(0, type = "n", xlim = limx, ylim = limy, asp = 1, frame.plot = FALSE, 
        axes = FALSE, ylab = ylab,  sub=sub, main = main,
        xlab = xlab )
    max.len <- max(strwidth(rownames(factors)), strwidth("abc"))/1.8
    limx = c(-max.len/2, x.max)
    cex <- min(cex, 20/x.max)
    for (v in 1:nvar) {
        var.rect[[v]] <- dia.rect(0, nvar - v + 1, rownames(factors)[v], 
            xlim = limx, ylim = limy, cex = cex, ...)
    }
    f.scale <- (nvar + 1)/(num.factors + 1)
    f.shift <- nvar/num.factors
    for (f in 1:num.factors) {
        fact.rect[[f]] <- dia.ellipse(0.6 * x.max, (num.factors + 
            1 - f) * f.scale, colnames(factors)[f], xlim = limx, 
            ylim = limy, e.size = e.size, ...)
  
#          text(x.max-1,(num.factors + 1 - f) * f.scale, 
#               labels=round(prop[f,2],digits) )
          text(x.max,(num.factors + 1 - f) * f.scale, 
              labels=round(prop[f,3],digits) )
            axis(2,at=nvar:1,labels=round(unique,digits),
                 las=1,tick=FALSE,font=1)

        for (v in 1:nvar) {
            if (simple && (abs(factors[v, f]) == max(abs(factors[v,]))) && 
                (abs(factors[v, f]) > cut) | 
                (!simple && (abs(factors[v, f]) > cut))) {
                if (fa) {
                dia.arrow(from = fact.rect[[f]], to = var.rect[[v]]$right, 
                  labels = round(factors[v, f], digits),
                   col = col[(sign(factors[v, f]) < 0) + 1], 
                   lty = ((sign(factors[v, f]) < 0) + 1))
                } else {
                dia.arrow(from = var.rect[[v]]$right, to = fact.rect[[f]], 
                  labels = round(factors[v, f], digits),
                   col = col[(sign(factors[v, f]) < 0) + 1], 
                   lty = ((sign(factors[v, f]) < 0) + 1))
                }
            }
        }
    }
    if (!is.null(Phi)) {
        for (i in 2:num.factors) {
            for (j in 1:(i - 1)) {
                if (abs(Phi[i, j]) > cut) {
                  dia.curve(from = fact.rect[[j]]$right, 
                            to = fact.rect[[i]]$right, 
                    labels = round(Phi[i, j], digits), 
                    scale = (i - j), ...)
               }
            }
        }
    }
    if (errors) {
        for (v in 1:nvar) {
            dia.self(location = var.rect[[v]], scale = 0.5, side = side)
        }
    }

  if (!is.null(addF)) {
     for (i in 1:nrow(addF)) {
       f <- addF[i,2]
       v <- addF[i,1]
       bond <- addF[i,3]
       lm <- addF[i,4]
       if(lm==0){
                dia.arrow(from = fact.rect[[f]], to = var.rect[[v]]$right, 
                  labels = round(bond, digits),
                   col = col[(sign(bond) < 0) + 1], 
                   lty = ((sign(bond) < 0) + 1))
       } else {
                dia.arrow2(from = fact.rect[[f]], to = var.rect[[v]]$right, 
                  labels = round(bond, digits),lmove=lm,
                   col = col[(sign(bond) < 0) + 1], 
                   lty = ((sign(bond) < 0) + 1))
       }
     }
  }
  if (!is.null(addP)) {
     for (i in 1:nrow(addP)) {
       f <- addP[i,2]
       v <- addP[i,1]
       bond <- addP[i,3]
       lm <- addP[i,4]
       if(lm==0){
                dia.arrow(from = var.rect[[v]]$right, to = fact.rect[[f]], 
                  labels = round(bond, digits),
                   col = col[(sign(bond) < 0) + 1], 
                   lty = ((sign(bond) < 0) + 1))
       } else {
                dia.arrow2(from = var.rect[[v]]$right, to = fact.rect[[f]], 
                  labels = round(bond, digits),lmove=lm,
                   col = col[(sign(bond) < 0) + 1], 
                   lty = ((sign(bond) < 0) + 1))
       }
     }
  }

}       ###################### end pca.diagram

if (is.null(protate)) {
  protate <- rotate
}
fres <- NULL
r.scores <- NULL
r.loadings <- NULL
r.eig <- NULL
r.names <- NULL

if (pm=="Facto" & is.null(n.obs))  {                     
  res <- PCA(data, ncp=ncp, graph=FALSE)
  load <- res$var$coord
  prop <- res$eig
  if (graph) {
    plot.PCA2(res,choix="var",atype=2,new.plot=new.plot,graph=TRUE)
  }

} else {
  res <- principal(data,nfactors=min(max.cp, ncol(data),ncp),n.obs=n.obs,
    rotate=protate,plot.scores=is.null(n.obs))
  load <- r.loadings <- res$loadings
  if (!is.null(res$Phi)) 
    Phi <- res$Phi
  prop <- colSums(load^2)
  prop <- cbind(prop,100*prop/ncol(data))
  prop <- r.eig <- cbind(prop, cumsum(prop[,2]))
  r.scores <- res$scores
  r.names <- names(res$R2)
}

pca.diagram(load, prop, Phi = Phi, sort = sort, labels = labels, 
    cut = cut, simple = simple, errors = errors, digits = digits, 
    e.size = e.size, rsize = rsize, side = side, main = main, cex = cex, 
    col2=col2,xlab=xlab, ylab=ylab, sub=sub, data.name=data.name, fa=FALSE,
    n.obs=n.obs,fm=fm, rotate=rotate, new.plot=new.plot, max.cp=max.cp, 
    addP=addP,addF=NULL, lmove=lmove, protate=protate,...)

if (fa) {
    fres <- fa(data,fm=fm,nfactors=nfactors, n.obs=n.obs,rotate=rotate)
 load <- fres$loadings
  if (!is.null(fres$Phi)) Phi <- fres$Phi
  prop <- colSums(load^2)
  prop <- cbind(prop,100*prop/ncol(data))
  prop <- cbind(prop, cumsum(prop[,2]))
 pca.diagram(fres, prop, fa=TRUE, Phi=Phi, sort = sort, labels = labels, 
    cut = cut, simple = simple, errors = TRUE, digits = digits, cex=cex,
    col2=col2, e.size = e.size, rsize = rsize, side = side, 
    main = main,data.name=data.name,protate=protate,
    sub=sub, ylab=ylab, xlab=xlab,
    fm=fm, rotate=rotate,n.obs=n.obs,
    new.plot=new.plot, max.cp=max.cp, 
    addF=addF, addP=NULL, lmove=lmove, ...)
}

return(invisible(list(res=res, fres=fres,r.names=r.names,
  r.eig=r.eig, r.loadings=r.loadings, r.scores=r.scores)))

}  # end PCA.diagram