/* GAUSS implementation of 
** Gradient projection for rotation in factor analysis
** 
** Translated into GAUSS
** from the original Matlab code written by
** 
** Bernaards and Jennrich downloadable from
** http://www.stat.ucla.edu/research/gpa/
** 
** For details, see Bernaards, C.A. and Jennrich, R.I. (in press) 
** Gradient Projection Algorithms and Software for Arbitrary 
** Rotation Criteria in Factor Analysis. 
** Educational and Psychological Measurement.
**
** and  
** 
** Jennrich, R.I. (2001). A simple general procedure for orthogonal 
** rotation. Psychometrika, 66, 289-306.
** Jennrich, R.I. (2002). A simple general method for oblique 
** rotation. Psychometrika, 67, 7-19.  
**
** 03/21/05 by Li Cai */;


/* These are commonly used utilities */;
proc Dg(x);
  retp(diagrv(eye(rows(x)),x));
endp;

proc tr(m);
  retp(sumc(diag(m)));
endp;

proc Jay(m);
  retp(ones(m,m));
endp;

proc sum(m);
  retp(sumc(vec(m)));
endp;

/* quartimin */;
proc(2) = vgQ_quartimin(L);
  local L2,k,q,M,Gq;
  L2 = L^2;
  k = cols(L);
  M = Jay(k)-eye(k);
  q = sum(L2.*(L2*M))/4;
  Gq = L.*(L2*M);
  retp(q,Gq);
endp;

/* Oblimin family */;
proc(2) = vgQ_oblimin(L,gam);
  local L2,k,p,M,q,Gq;
  /*
  if (gam == 0) "Oblimin Quartimin"
  if (gam == .5) "Oblimin Biquartimin"
  if (gam == 1) "Oblimin Covarimin"
  */;
  L2 = L^2;
  k = cols(L);
  p = rows(L);
  M = Jay(k)-eye(k);

  q = sum(L2.*((eye(p)-gam.*(1/p).*Jay(p))*L2*M))/4;
  Gq = L.*((eye(p)-gam.*(1/p).*Jay(p))*L2*M);
  retp(q,Gq);
endp;

/* Target */;
proc(2) = vgQ_target(L,Target);
  local q,Gq;
  q = sum((L-Target)^2);
  Gq = 2.*(L-Target);
  retp(q,Gq);
endp;

/* Partially specified target */;
proc(2) = vgQ_pst(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) */;
  local Btilde,q,Gq;
  Btilde = W.*Target;
  q = sum((W.*L-Btilde)^2);
  Gq = 2.*(W.*L-Btilde);
  retp(q,Gq);
endp;

/* Oblimax */;
proc(2) = vgQ_oblimax(L);
  local q,Gq,L4,L3,L2;
  L4 = L^4;
  L3 = L^3;
  L2 = L^2;
  q = -(ln(sum(L4))-2*ln(sum(L2)));
  Gq = -(4.*L3./sum(L^4)-4.*L./sum(L2));
  retp(q,Gq);
endp;

/* Minimum entropy */;
proc(2) = vgQ_entropy(L);
  local q,Gq,L2,W;
  L2 = L^2;
  W = L2 .eq 0;
  q = -sum(L2.*ln(L2+W))/2;
  Gq = -(L.*ln(L2+W)+L);
  retp(q,Gq);
endp;

/* Quartimax */;
proc(2) = vgQ_quartimax(L);
  local L2,q,Gq;
  L2 = L^2;
  q = -sum(diag(L2'*L2))/4;
  Gq = -L^3;
  retp(q,Gq);
endp;

/* Varimax */;
proc(2) = vgQ_varimax(L);
  local QL,means,L2,q,Gq;
  L2 = L^2;
  means = meanc(L2)';
  QL = L2-means;
  q = -sqrt(sum(diag(QL'*QL)))^2/4;
  Gq = -L.*QL;
  retp(q,Gq);
endp;

/* Simplimax */;
proc(2) = vgQ_simplimax(L,k);
  /* m: Number of close to zero loadings 
     if k = 0 then m = rows(L) */;
  local Imat,L2,q,Gq,m,sorted;
  m = k;
  if (k eq 0);
    m = rows(L);
  endif;
  L2 = L^2;
  sorted=sortc(vec(L2),1);
  Imat = L^2 .le sorted[m];
  q = sum(Imat.*L2);
  Gq = 2.*Imat.*L;
  retp(q,Gq);
endp;

/* Geomin */;
_geomin_eps = .01;
proc(2) = vgQ_geomin(L);
  local k,p,L2,pro,q,Gq,eps;
  eps = _geomin_eps;
  k = cols(L);
  p = rows(L);
  L2 = L^2+eps;
  pro = exp((sumc(ln(L2')))/k);
  q = sumc(pro);
  Gq = (2/k).*(L./L2).*(pro.*ones(p,k));
  retp(q,Gq);
endp;

/* Crawford-Ferguson family */;
proc(2) = vgQ_cf(L,kappa);
/*   
  kappa <- 0 # Quartimax 
  kappa <- 1/p # Varimax
  kappa <- m/(2*p) # Equamax
  kappa <- (m-1)/(p+m-2) # Parsimax
  kappa <- 1 # Factor parsimony
*/;
  local k,p,N,M,L2,f1,f2,q,Gq;
  k = cols(L);
  p = rows(L);
  N = Jay(k)-eye(k);
  M = Jay(p)-eye(p);
  L2 = L^2;
  f1 = (1-kappa)*sumc(diag(L2'*(L2*N)))/4;
  f2 = kappa*sumc(diag(L2'*(M*L2)))/4;
  q = f1 + f2;
  Gq = (1-kappa).*L.*(L2*N) + kappa.*L.*(M*L2);
  retp(q,Gq);
endp;

proc(4) = GPForth(A,Ts,MAXCC,MAXI,Tol);
  local al,L,f,ft,Gq,G,iter,EC,M,S,Gp,norm,i,X,U,D,V,Tt,T,Lh,Th;
  EC = MAXCC;
  al=1;
  L=A*Ts;
  {f,Gq} = vgQ_quartimax(L);
  G = A'*Gq;
  T = Ts;
  for iter (1,MAXCC,1);
    M=T'*G;
    S=(M+M')/2;
    Gp=G-T*S;
    norm=sqrt(sum(Gp^2));
    if (norm lt Tol);
      EC = iter;
      break;
    endif;
    al = 2*al;
    for i (1,MAXI,1);
      X=T-al*Gp;
      {U,D,V}=svd1(X);
      Tt=U*V';
      L=A*Tt;
      {ft,Gq}=vgQ_quartimax(L);
      if (ft lt (f-.5*(norm^2)*al));
        break;
      endif;
      al=al/2;
    endfor;
    T=Tt;
    f=ft;
    G=A'*Gq;
  endfor;
  Th=T;
  Lh=A*T;
  retp(Th,Lh,f,EC); 
endp;

/* test out the orthogonal routine using Harman's example */;
A = {.830 -.396,
     .818 -.469,
     .777 -.470,
     .798 -.401,
     .786  .500,
     .672  .458,
     .594  .444,
     .647  .333 };

{Th,Lh,f,EC} = GPForth(A,eye(2),500,10,1e-5);

proc(5) = GPForb(A,Ts,MAXCC,MAXI,Tol);
  local al,T,Ti,L,f,Gq,G,iter,Gp,norm,i,X,v,Tt,ft,Th,Lh,Phi,EC;
  EC = MAXCC;
  al=1;
  T=Ts;
  Ti=inv(T);
  L=A*Ti';
  {f,Gq}=vgQ_quartimin(L);
  G=-(L'*Gq*Ti)';
  for iter (1,MAXCC,1);
    Gp=G-T*Dg(sumc(T.*G));
    norm=sqrt(sum(Gp^2));
    if (norm lt Tol);
      EC = iter;
      break;
    endif;
    al=2*al;
    for i (1,MAXI,1);
       X=T-al*Gp;
       v=1./sqrt(sumc(X^2));
       Tt=X*Dg(v);
       Ti=inv(Tt);
       L=A*Ti';
       {ft,Gq}=vgQ_quartimin(L);
       if (ft lt (f-.5*(norm^2)*al));
         break;
       endif;
       al=al/2;
    endfor;
    T=Tt;
    f=ft;
    G=-(L'*Gq*Ti)';
  endfor;
  Th=T;
  Lh=L;
  Phi=T'*T;
  retp(Th,Lh,Phi,f,EC);
endp;

/* test out the oblique routine using Harman's example */;
A = {.830 -.396,
     .818 -.469,
     .777 -.470,
     .798 -.401,
     .786  .500,
     .672  .458,
     .594  .444,
     .647  .333 };

{Th,Lh,Phi,f,EC}=GPForb(A,eye(2),500,10,1e-5);

/* Generate random start matrix */;
proc RandomStart(k);
  local m,q,r;
  m = rndn(k,k);
  {q,r} = qqr(m);
  retp(q);
endp;