# This file contains editable Splus/R code for general GP
# algorithms for rotation. See "http://www.stat.ucla.edu/research"
# for a discussion of these.

# GPForth.df is the main GP algorithm for orthogonal rotation.
# GPFoblq.df is the main GP algorithm for oblique rotation.
# Gf computes the numerical derivative and is needed by GPForth.df and GPFoblq.df.
# For both algorithms is required: a loadings matrix A. Optional  
# a initial rotation matrix Tmat. By default this is the identity matrix.
# Optional: the rotation method to be used. Between quation marks have to
# be the last part of the name of the ff function, e.g. for ff.varimax
# the argument is "varimax". Identical arguments can be used for oblique
# rotation. Some rotation criteria (including simplimax, pst, procrustes, 
# cf,...) require one or more additional arguments. For example, simplimax
# needs the number of 'close to zero loadings'. This is given included as
# the extra argument k=27. Check out the rotation methods for details. 
# When a new rotation method is implemented, and it needs an additional
# argument then this is the easiest way to pass it to the function. 
#
# New rotation methods need to be programmed as ff.newmethod. The only
# inputs are the matrices Tmat and A, and potential additional arguments. The
# output consists of the value f  of the criterion.

GPForth.df <- function(A,Tmat=diag(ncol(A)),method="varimax",...){
  Gf <- function(Tmat, A, method,...){
   k <- nrow(Tmat)
   ep <- .0001
   G <- Z <- matrix(0,k,k)
   for (r in 1:k){
     for (s in 1:k){
        dT <- Z
        dT[r,s] <- ep
        p1 <- get(paste("ff",method,sep="."))(A %*% (Tmat+dT),...)
        p2 <- get(paste("ff",method,sep="."))(A %*% (Tmat-dT),...)
        G[r,s] <- (p1 - p2)/(2*ep)
      }
    }
  return(G)
  }
 al <- 1
 Table <- NULL
 L <- A %*% Tmat
 for (iter in 0:100){
   f <- get(paste("ff",method,sep="."))(L,...)
   G <- Gf(Tmat, A, method, ...)
   M <- crossprod(Tmat,G)
   S <- (M + t(M))/2
   Gp <- G-Tmat %*% S
      s <- sqrt(sum(diag(crossprod(Gp))))
   Table <- rbind(Table,c(iter,f,log10(s),al))
   if (s < 1e-5)
         break
   al <- 2*al
   for (i in 0:10){
     X <- Tmat - al * Gp
     UDV <- svd(X)
     Tmatt <- UDV$u %*% t(UDV$v)
     L <- A %*% Tmatt
     ft <- get(paste("ff",method,sep="."))(L,...)
     if (ft < (f-.5*s^2*al))
       break
     al <- al/2
    }
    Tmat <- Tmatt
  }
  Th <- Tmat
  Lh <- A %*% Tmat
  return(list(Lh=Lh,Th=Th,Table=Table,method=method,orthogonal=T))
}
 


GPFoblq.df <- function(A,Tmat=diag(ncol(A)),method="quartimin",...){
  Gf <- function(Tmat, A, method,...){
   k <- nrow(Tmat)
   ep <- .0001
   G <- Z <- matrix(0,k,k)
   for (r in 1:k){
     for (s in 1:k){
        dT <- Z
        dT[r,s] <- ep
        p1 <- get(paste("ff",method,sep="."))(A %*% t(solve(Tmat+dT)),...)
        p2 <- get(paste("ff",method,sep="."))(A %*% t(solve(Tmat-dT)),...)
        G[r,s] <- (p1 - p2)/(2*ep)
      }
    }
  return(G)
  }
 al <- 1
 Table <- NULL
 L <- A %*% t(solve(Tmat))
 for (iter in 0:500){
    f <- get(paste("ff",method,sep="."))(L,...)
    G <- Gf(Tmat, A, method)
    Gp <- G-Tmat %*% diag(apply(Tmat*G,2,sum))
    s <- sqrt(sum(diag(crossprod(Gp))))
    Table <- rbind(Table,c(iter,f,log10(s),al))
    if (s < 1e-5)
       break
     al <- 2*al
     for (i in 0:10){
       X <- Tmat-al*Gp
       v <- 1/sqrt(apply(X^2,2,sum))
       Tmatt <- X %*% diag(v)
       L <- A %*% t(solve(Tmatt))
       ft <- get(paste("ff",method,sep="."))(L,...)
       if (ft < (f-.5*s^2*al))
          break
       al <- al/2
     } 
     Tmat <- Tmatt
   }
   Th <- Tmat
   Lh <- A %*% t(solve(Tmat))
   Phi <- t(Tmat) %*% Tmat
   return(list(Lh=Lh,Phi=Phi,Th=Th,Table=Table,method=method,orthogonal=F))
}


ff.quartimax <- function(L){
  L2 <- L * L
  f <- -sum(L2 * L2)/4
  return(f)
}

ff.varimax <- function(L){
  L2 <- L * L
  QL <- sweep(L^2,2,apply(L^2,2,mean),"-")
  f <- -sqrt(sum(diag(crossprod(QL))))^2/4
  return(f)
}

ff.quartimin <- function(L){
  L2 <- L^2
  k <- ncol(L)
  M <- matrix(1,k,k)-diag(k)
  f <- sum(L2 * (L2 %*% M))/4
  return(f)
} 

ff.oblimin <- function(L,gam=0){
  # if (gam == 0) Method <- "Oblimin Quartimin"
  # if (gam == .5) Method <- "Oblimin Bi-quartimin"
  # if (gam == 1) Method <- "Oblimin Covarimin"
  k <- ncol(L)
  p <- nrow(L)
  N <- matrix(1,k,k)-diag(k)
  f <- sum(L^2 * (diag(p)-gam*matrix(1/p,p,p)) %*% L^2 %*% N)/4
  return(f)
}

ff.target <- function(L,Target){
  # Needs Target matrix, e.g.  Target <-matrix(c(rep(9,4),rep(0,8),rep(9,4)),8) 
  f <- sum((L-Target)^2)
  return(f)
}

ff.pst <- function(L,W,Target){
  # Needs weight matrix W with 1's at specified values, 0 otherwise
  # e.g. W = matrix(c(rep(1,4),rep(0,8),rep(1,4)),8). 
  # When W has only 1's this is procrustes rotation
  # Needs a Target matrix Target with hypothesized factor loadings.
  # e.g. Target = matrix(0,8,2)
  Btilde <- W * Target
  f <- sum((W*L-Btilde)^2)
  return(f)
}

ff.oblimax <- function(L){
  f <- -(log(sum(L^4))-2*log(sum(L^2)))
  return(f)
}

ff.entropy <- function(L){
  f <- -sum(L^2 * log(L^2 + (L^2==0)))/2
  return(f)
}
 
ff.cubimax <- function(L){
  f <- -sum(diag(t(L^2) %*% abs(L)))
  return(f)
}

ff.simplimax <- function(L,k=nrow(L)){
  # k: Number of close to zero loadings
  Imat <- sign(L^2 <= sort(L^2)[k])
  f <- sum(Imat*L^2)
  return(f)
}

ff.ss <- function(L){
  # Method <- "Forced Simple Structure"
  # m: Number of close to zero loadings
  m <- ncol(L)
  zm <- m+2
  Imat <- matrix(0,nrow(L),ncol(L))
  for (i in 1:nrow(L))
    Imat[i,abs(L[i,]) <= sort(abs(L[i,]))[1]] <- 1
  for (j in 1:m)
    Imat[abs(L[,j]) <= sort(abs(L[,j]))[zm],j] <- 1
  for (i in 1:(m-1)){
    for (j in (i+1):m){
      nz <- sum(Imat[,i] + Imat[,j]==1)
      while (nz < zm){
	tbc <- c(abs(L[,i]), abs(L[,j]))
        lo <- sum(c(Imat[,i], Imat[,j])==1)
	tbc[sort(abs(tbc))[lo+1]] <- 1
	nz <- sum(Imat[,i] + Imat[,j]==1)
      }
    }
  }
  f <- sum(Imat*L^2)
  return(f)
}

ff.bentler <- function(L){
  L2 <- L^2
  M <- crossprod(L2)
  D <- diag(diag(M))
  f <- -(log(det(M))-log(det(D)))/4
  return(f)
}

ff.tandemI <- function(L){  # Tandem Criterion I, Comrey, 1967.
  LL <- (L %*% t(L))
  LL2 <- LL^2
  f <- -sum(diag(crossprod(L^2, LL2 %*% L^2)))
  return(f)
}

ff.tandemII <- function(L){  # Tandem Criterion II, Comrey, 1967.
  LL <- (L %*% t(L))
  LL2 <- LL^2
  f <- sum(diag(crossprod(L^2, (1-LL2) %*% L^2)))
  return(f)
}

ff.geomin <- function(L,eps=.01){
  k <- ncol(L)
  L2 <- L^2+eps
  pro <- exp(apply(log(L2),1,sum)/k) #apply(L2,1,prod)^(1/k)
  f <- sum(pro)
  return(f)
}

ff.infomax <- function(L){
  k <- ncol(L)
  S <- L^2
  s <- sum(S)
  s1 <- apply(S, 1, sum)
  s2 <- apply(S, 2, sum)
  E <- S/s
  e1 <- s1/s
  e2 <- s2/s
  Q0 <- sum(-E * log(E))
  Q1 <- sum(-e1 * log(e1))
  Q2 <- sum(-e2 * log(e2))
  f <- log(k) + Q0 - Q1 - Q2
  return(f)
}

ff.mccammon <- function(L){
  S <- L^2
  s <- sum(S)
  s2 <- apply(S, 2, sum)
  e0 <- sweep(S, 2, s2, "/")
  e2 <- s2/s
  Q0 <- sum(e0 * log(e0))
  Q2 <- sum(e2 * log(e2))
  f <- Q0/Q2
  return(f)
}

ff.cf <- function(L,kappa=0){
  k <- ncol(L)
  p <- nrow(L)
  # kappa <- 0 # Quartimax 
  # kappa <- 1/p # Varimax
  # kappa <- m/(2*p) # Equamax
  # kappa <- (m-1)/(p+m-2) # Parsimax
  # kappa <- 1 # Factor parsimony
  # Method <- paste("Crawford-Ferguson:k=",kappa,sep="")
  N <- matrix(1,k,k)-diag(k)
  M <- matrix(1,p,p)-diag(p)
  L2 <- L^2
  f1 <- (1-kappa)*sum(diag(crossprod(L2,L2 %*% N)))/4
  f2 <- kappa*sum(diag(crossprod(L2,M %*% L2)))/4
  f <- f1 + f2
  return(f)
}

##
## UTILITY FUNCTIONS
##

# Kaisers standardized variables
# Can be used for standardization and de-standardization.
# For standardization: only provide 1 argument: the matrix A
# output: standardized A and weights W
# For destandardization: provide 2 arguments: matrix A
# and weights W.
# output: destandardized A and weights W.
Standardize <- function(A,W){
 if (nargs()==1)
   W <- sqrt(apply(A^2,1,sum))
 else 
   W <- 1/W
 A <- sweep(A,1,W,"/")
 return(A,W)
}

# Random Start
# k: number of dimensions
# orthogonal random start? (Yes, by default).
Random.Start <- function(k,orthogonal=T){
 mat <- matrix(rnorm(k*k),k)
 if (orthogonal){
   ans <- qr.Q(qr(mat))
 }
 else{
   ans <- mat %*% diag(sqrt(1/diag(t(mat) %*% mat)))
 }
 ans
}