/*
     GPFOBL: oblique rotation for factor analysis
     by: HervŽ M. CACI (hcaci@wanadoo.fr)

     version 1.0 -- Date: February 19th, 2003
     Translated from an SPSS macro for OBLIMIN only
     by Coen A. Bernaards & Robert I. Jennrich
     http://www.stat.ucla.edu/research

     version 1.1 -- Date: March 15th, 2003
     added GAMMA and FACTOR options

     version 1.2 -- Date: March 23rd, 2003
     bugs correction & other rotation criteria added

     version 1.3 -- Date: May 2nd, 2003
     added Hyperplane Count on rotated loadings
*/

program define gpfobl
version 5.0
*! Version 1.3 -- May 2nd, 2003

/* Parse standard STATA commands */
parse "`*'", parse(" ,")
confirm exist `1'
matrix A0=`1'
mac shift 1

/* Parse options */

#delimit ;
local options "Factor(int -1) OBLIMIN(real -2) OBLIMAX CFerguson(real -2) Geomin(real -2) IPS";
#delimit cr
parse "`*'"

/* Extensive check of the rotation criterion */

di in gr _n "Gradient Projection Algorithm for arbitrary rotation criteria" _n
local subr=0
if "`oblimax'" != "" { di in gr _n "OBLIMAX criterion" _n _dup(17) "-"
                       local subr=1 }
if `oblimin' != -2 {
  if `subr'!=0 { di in red "More than one criterion specified"
                 exit }
  if `oblimin' == .50 { di in gr _n "OBLIMIN rotation" _n _dup(16) "-"
                        local gamma = 0.50 }
    else if `oblimin' == 0 { di in gr _n "QUARTIMIN rotation" _n _dup(18) "-" _n
                             local gamma=0.00 }
    else if `oblimin' == 1 { di in gr _n "COVARIMIN rotation" _n _dup(18) "-" _n
                             local gamma = 1 }
    else if `oblimin'<-1.00 | `oblimin'>1.00 { di in red "Gamma is out of range."
                                               exit }
    else { di in gr "OBLIMIN rotation" _n _dup(16) "-"
           local gamma=`oblimin'
           di in gr "Gamma parameter = " in ye %5.4f `gamma' 
  }
  local subr=2
}
if `geomin' != -2 {
  if `subr'!=0 { di in red "More than one criterion specified."
                 exit
  }
  di in gr _n "GEOMIN criterion" _n _dup(16) "-"
  if `geomin'<-1.00 | `geomin'>1.00 { di in red "Epsilon is out of range."
                                      exit }
  local subr=3
  local epsilon=`geomin'
  di in gr "Epsilon parameter = " in ye %5.4f `epsilon'
}
if "`ips'"!="" { di in gr _n "Bentler's Invariant Pattern Simplicity"
                 di in gr _dup(38) "-"
                 if `subr'!=0 { di in red "More than one criterion specified."
                                exit }
                 local subr=4
}

if `cferguson' != -2 {
  if `subr'!=0 { di in red "More than one criterion specified."
                 exit
  }
  di in gr _n "Crawford-Ferguson criterion" _n _dup(27) "-"
  if `cferguson'==0 { di in gr _n "Kappa = 0 -> Quartimax" }
  if `cferguson'==1 { di in gr _n "Kappa = 1 -> Factor parsimony" }
  if `cferguson'<-1.00 | `cferguson'>1.00 { in red "Kappa is out of range."
                                            exit }
    else { di in gr "Kappa = " in ye %3.2f `cferguson' }
  local kappa=`cferguson'
  local subr=5
}

if `subr'==0 { di in red "No criterion specified."
               exit }

/* To date, the Tmat matrix is IDENTITY... may change in the future */

di in gr "(Note: Initial normal rotation matrix is IDENTITY by default)."
matrix Tmat=I(colsof(A0))

/* FACTOR option */

if `factor'==-1 { local factor=colsof(A0) }
if `factor'>colsof(A0) {
    di in red "Too many factors specified"
    exit
  }
if `factor'<2 {
    di in red "Not enough factors specified"
    exit
  }
matrix A=A0[.,1..`factor']
matrix Tmat=Tmat[1..`factor',1..`factor']

/* Main program starts here */

local maxiter=500
local a1=1
mat Ti=inv(Tmat)
mat L=A*Ti'
local k=colsof(L)         /* Number of factors   */
local p=rowsof(L)         /* Number of variables */
di in gr _n "Number of Factors  : " %3.0f in ye `k'
di in gr "Number of Variables: " %3.0f in ye `p'
if `k'<2 {di in red _n "Too small number of factors..."
          exit}
if `k'>=`p' {di in red _n "Not enough variables for that number of factor..."
             exit}

/*
  Call the appropriate subroutine
*/

global ft=0
mat Gq=J(`p',`k',0)
if `subr'==1 { oblimax L Gq }
 else if `subr'==2 { oblimin L `gamma' Gq }
   else if `subr'==3 { geomin L `epsilon' Gq }
     else if `subr'==4 { ips L Gq }
       else if `subr'==5 { cferg L `kappa' Gq }

/*
  Long loop
*/

local f=$ft
matrix G=-(L'*Gq*Ti)'
local iter=0
while `iter'<=`maxiter' {
  matrix csum=J(1,`k',0)
  matrix TG=J(`k',`k',0)
  local i=1
  while `i'<=`k' {
    local j=1
    while `j'<=`k' {
      matrix TG[`i',`j']=Tmat[`i',`j']*G[`i',`j']
      matrix csum[1,`j']=csum[1,`j']+TG[`i',`j']
      local j=`j'+1
    }
    local i=`i'+1
  }
  matrix Gp=G-Tmat*diag(csum)
  matrix Y=Gp'*Gp
  local s=sqrt(trace(Y))
/*  di "`iter' | $ft | `s1' | `a1'" */
  if int(`iter'/5)==(`iter'/5) { di in ye "." _c }
  local s1=log10(`s')
  if `s'<0.00001 {
     di in gr _n "*** Algorithm converged after " %3.0f in ye `iter' /*
     */ in gr " iterations ***"
     local iter=`maxiter'+10
  }
  else {
    local a1=`a1'*2
    local i=0
    while `i'<11 {
    matrix X=J(`k',`k',0)
    matrix X2=J(`k',`k',0)
    matrix csum=J(1,`k',0)
      local x=1     /* This loop computes X=Tmat-a1&*Gp */
      while `x'<=`k' {
        local y=1
        while `y'<=`k' {
          matrix X[`x',`y']=Tmat[`x',`y']-(`a1'*Gp[`x',`y'])
          matrix X2[`x',`y']=X[`x',`y']*X[`x',`y']
          matrix csum[1,`y']=csum[1,`y']+X2[`x',`y']
          local y=`y'+1
        }
        local x=`x'+1
      }
      mat v=J(1,`k',0)
      local y=1
      while `y'<=`k' {
        matrix v[1,`y']=1/sqrt(csum[1,`y'])
        local y=`y'+1
      }
      matrix Tt=X*diag(v)
      matrix Ti=inv(Tt)
      matrix L=A*Ti'

/*
   Call the appropriate subroutine
*/

      global ft=0
      mat Gq=J(`p',`k',0)
      if `subr'==1 { oblimax L Gq }
        else if `subr'==2 { oblimin L `gamma' Gq }
          else if `subr'==3 { geomin L `epsilon' Gq }
            else if `subr'==4 { ips L Gq }
              else if `subr'==5 { cferg L `kappa' Gq }

/* Common part continues... */

      matrix G=-(L'*Gq*Ti)'
      local f0=`f'-(`s'*`s'*`a1')/2
      if $ft>`f0' { local a1=`a1'/2 }
        else if $ft<`f0' { local i=10 } /* Break the loop... */
      local i=`i'+1
    } /* End of the `i' loop */
  matrix Tmat=Tt
  local f=$ft
  matrix G=-(L'*Gq*Ti)'
  } /* end of "else" statement... */
  local iter=`iter'+1
}
if `iter'==`maxiter'+1 {
  di in red _n "*** Algorithm didn't converge after `maxiter' iterations." }
else {
  di in gr _n "Transformation matrix"
  mat list Tmat, format(%4.3f) noheader
  di in gr _n "Factor loadings after rotation"
  mat list L, format(%4.3f) noheader
  di in gr _n "Factor intercorrelations"
  mat phi=Tmat'*Tmat
  mat list phi, format(%4.3f) noheader
}

/* Compute the hyperplane counts per factor */

mat HCOUNT=J(1,`k',0)
local i=1
while `i'<=`k' {
  local j=1
  local s=0
  while `j'<=`p' {
    if abs(L[`j',`i'])<=.100 { local s=`s'+1 }
    local j=`j'+1
  }
  mat HCOUNT[1,`i']=`s'
  local i=`i'+1
}
di in gr _n "Hyperplane count: loadings in range [-0.10 ; +0.10]"
mat list HCOUNT,format(%3.0f) noheader

end

/*
   OBLIMIN subroutine: Matrixname gamma Gq
*/

program define oblimin
local k=colsof(`1')         /* Number of factors   */
local p=rowsof(`1')         /* Number of variables */
mat N=J(`k',`k',1)-I(`k')
mat L2=J(`p',`k',0)
mat Y=J(`p',`p',-`2'/`p')
local x=1
while `x'<=`p' {
  mat Y[`x',`x']=1-(`2'/`p')
  local y=1
  while `y'<=`k' {
    mat L2[`x',`y']=L[`x',`y']*L[`x',`y']
    local y=`y'+1
  }
  local x=`x'+1
}
mat Y=Y*L2*N
mat F=J(`p',`k',0)
mat G=J(`p',`k',0)
local x=1
local f1=0
while `x'<=`p' {
  local y=1
  while `y'<=`k' {
    mat F[`x',`y']=L2[`x',`y']*Y[`x',`y']
    local f1=`f1'+F[`x',`y']
    mat G[`x',`y']=L[`x',`y']*Y[`x',`y']
    local y=`y'+1
  }
  local x=`x'+1
}
global ft=`f1'/4
matrix `3'=G
end

/*
   OBLIMAX subroutine: Matrixname Gq
*/

program define oblimax
local k=colsof(`1')         /* Number of factors   */
local p=rowsof(`1')         /* Number of variables */
local i=1
local L2=0         /* msum(L&*L)  */
local L4=0         /* msum(L&**4) */
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    local L2=`L2'+`1'[`i',`j']^2
    local L4=`L4'+`1'[`i',`j']^4
    local j=`j'+1
  }
  local i=`i'+1
}
global ft=-(log(`L4')-2*log(`L2'))
mat G=J(`p',`k',0)
local i=1
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    mat G[`i',`j']=-(4*(`1'[`i',`j']^3)/`L4'-4*`1'[`i',`j']/`L2')
    local j=`j'+1
  }
  local i=`i'+1
}
matrix `2'=G
end

/*
   GEOMIN subroutine: Matrixname epsilon Gq
*/

program define geomin
local k=colsof(`1')         /* Number of factors   */
local p=rowsof(`1')         /* Number of variables */
local f1=0
mat L2=J(`p',`k',0)
mat LL2=J(`p',`k',0)
mat SL=J(`p',1,0)
local i=1
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    mat L2[`i',`j']=(`1'[`i',`j']*`1'[`i',`j'])+`2'
    mat LL2[`i',`j']=log(L2[`i',`j'])
    mat SL[`i',1]=SL[`i',1]+LL2[`i',`j']        /* RSUM(LL2) */
    local j=`j'+1
  }
  local f1=`f1'+exp(SL[`i',1]/`k')
  local i=`i'+1
}
global ft=`f1'
mat RM=LL2*J(`k',`k',1/`k')
mat G=J(`p',`k',0)
local i=1
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    mat RM[`i',`j']=exp(RM[`i',`j'])
    mat G[`i',`j']=(2*`k')*(`1'[`i',`j']/L2[`i',`j'])*RM[`i',`j']
    local j=`j'+1
  }
  local i=`i'+1
}
matrix `3'=G
end

/*
   Bentler's Invariant Pattern Simplicity
*/

program define ips
local k=colsof(`1')         /* Number of factors   */
local p=rowsof(`1')         /* Number of variables */
local f1=0
mat L2=J(`p',`k',0)
local i=1
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    mat L2[`i',`j']=`1'[`i',`j']*`1'[`i',`j']
    local j=`j'+1
  }
  local i=`i'+1
}
mat M=L2'*L2
mat D=diag(vecdiag(M))
global ft=-(log(det(M))-log(det(D)))/4
mat G=J(`p',`k',0)
mat X=L2*(inv(M)-inv(D))
local i=1
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    mat G[`i',`j']=-`1'[`i',`j']*X[`i',`j']
    local j=`j'+1
  }
  local i=`i'+1
}
mat `2'=G
end

/*
   CRAWFORD-FERGUSON: kappa = 0             -> Quartimax
                      kappa = 1/p           -> Varimax
                      kappa = k/(2p)        -> Equamax
                      kappa = (k-1)/(p+k-2) -> Parsimax
                      kappa = 1             -> Factor parsimony
*/

program define cferg
local k=colsof(`1')         /* Number of factors   */
local p=rowsof(`1')         /* Number of variables */
local f1=0
mat L2=J(`p',`k',0)
local i=1
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    mat L2[`i',`j']=`1'[`i',`j']*`1'[`i',`j']
    local j=`j'+1
  }
  local i=`i'+1
}
mat N=L2*(J(`k',`k',1)-I(`k'))
mat M=(J(`p',`p',1)-I(`p'))*L2
mat NN=L2'*N
mat MM=L2'*M
local f1=(1-`2')*trace(NN)/4
local f2=`2'*trace(MM)/4
global ft=`f1'+`f2'
mat G=J(`p',`k',0)
local i=1
while `i'<=`p' {
  local j=1
  while `j'<=`k' {
    mat G[`i',`j']=(1-`2')*`1'[`i',`j']*N[`i',`j'] + /*
*/ `2'*`1'[`i',`j']*M[`i',`j']
    local j=`j'+1
  }
  local i=`i'+1
}
mat `3'=G
end