%OBLIMIN        Direct Oblimin rotation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  INPUT
%               PP   ->  Loading matrix to be rotated
%               gg   ->  Gamma parameter
%               conv ->  Convergence value
%               nrs  ->  Number of random starts
% OUTPUT
%               Pr   ->  Rotated Loading matrix
%               TTr  ->  Inter-factors correlation matrix
%               OCR  ->  Oblimin criterium
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Pr, TTr, OCRr ] = oblimin (PP, gg,conv,nrs);

for nos=1:nrs

      [NV, NR] = size(PP);
      A      = PP*orth(rand(NR,NR) - .5);
      TOLCOS = 1.0-conv;
      GA     = gg/NV;
      S = diag(A*A');
      FL = zeros(NV,1);
      TT = eye(NR);
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SIMPLICITY CRITERIUM (F0)

      C = sum(A .* A);
      V = sum(A .* A .* A .* A);
      for J=1:NR,
        V(J)=V(J)-GA*C(J)*C(J);
      end;
      D = sum(C');
      G = sum(V');
      H = sum(S .* S);
      H=H-GA*D*D ;
      F0=H-G;
      
      L = 0; NROT = 1;
      while (NROT > 0),                            % LINE 60 EN BROWNE

        COMIN = 2;
        NROT  = 0;

        for IP=1:NR,
          for IQ=1:NR,
              if (IP ~= IQ),

                 DP=D; GP=G;
                 D=D-C(IP)-C(IQ);
                 G=G-V(IP)-V(IQ);
                 P=0; R=0; T=0; O=0; Y=0; Z=0;
                 for I=1:NV,
                   A1=A(I,IP);
                   A2=A(I,IQ);
                   AA=A1*A1;
                   BB=A2*A2;
                   AB=A1*A2;
                   FL(I)=S(I);
                   S(I)=S(I)-AA-BB;
                   Z=Z+AB;
                   R=R+AA*BB;
                   P=P+AB*AA;
                   T=T+AA*S(I);
                   O=O+AB*S(I);
                 end;                 
                 X=TT(IP,IQ);
                 GAC=GA*C(IP);
                 R=R-GAC*C(IQ);
                 P=P-GAC*Z;
                 O=O-GA*Z*D;
                 T=T-GAC*D;                 
                 P1=1.5*(X-P/V(IP));
                 Q1=0.5*(V(IP)-4*X*P+R+2*T)/V(IP);
                 R1=0.5*(X*(T+R)-P-O)/V(IP);                 
                 A2 = cubic(P1,Q1,R1,V(IP));
                 A22=A2*A2;
                 P=1+2*X*A2+A22;
                 if (P < 0), P = (-1) * P; end;
                 A1=sqrt(P);
                 A3=A2/A1;
                 ABCOS=abs((1+A2*X)/A1);
                 COMIN=min([COMIN ABCOS]);
                 if (ABCOS > TOLCOS),
                    D=DP; G=GP;
                    S = FL;
                 else                              % LINE 181 EN BROWNE
                    NROT=NROT+1;
                    A11=P*P;
                    C(IP)=P*C(IP);
                    V(IP)=A11*V(IP);
                    Z=0; Y=0;
                    for I=1:NV,
                      A(I,IQ)=A(I,IQ)-A2*A(I,IP);
                      A(I,IP)=A1*A(I,IP);
                      BB=A(I,IQ)*A(I,IQ);
                      Z=Z+BB*BB;
                      Y=Y+BB;
                      S(I)=S(I)+A(I,IP)*A(I,IP)+BB;
                    end;
                    V(IQ)=Z-GA*Y*Y;
                    C(IQ)=Y;
                    D=D+Y+C(IP);
                    G=G+V(IP)+V(IQ);
                    A1=1/A1;                   
                    for I=1:NR,
                       TT(I,IP)=A1*TT(I,IP)+A3*TT(I,IQ);
                       TT(IP,I)=TT(I,IP);
                    end;
                    TT(IP,IP)=1;
                 end;                              % LINE 81 EN BROWNE
              end;                                 % LINE 81 EN BROWNE
          end;                                     % LINE 81 EN BROWNE
        end;                                       % LINE 81 EN BROWNE
        
        C(NR) = sum(A(:,NR) .* A(:,NR));
        V(NR) = sum(A(:,NR) .* A(:,NR) .* A(:,NR) .* A(:,NR));
        V(NR) = V(NR) - GA * C(NR) * C(NR);

        D = sum(C');
        G = sum(V');
        S = diag(A*A');
        H = sum(S .* S);
        F = H - GA*D*D - G;
        if (NROT > 0), F0 = F; end;
      end;                                         % LINE 80 EN BROWNE
     
      for J=1:NR,
        C(J)=0.0;
        if (sum(A(:,J)) < 0),
           for I=1:NV,
               A(I,J) = (-1) * A(I,J);
           end;
           C(J) = 1;
        end;
      end;

      for I=1:NR,
        if (C(I) > 0),
           for J=1:NR,
              TT(I,J)= (-1) * TT(I,J);
              TT(J,I)= (-1) * TT(J,I);
           end;
        end;
      end;

  if nos == 1
    Fr   = F;  
    OCRr = F;
    Pr   = A;
    TTr  = TT;
  elseif F < Fr
    Fr   = F;
    OCRr = F;
    Pr   = A;
    TTr  = TT;
  end;

end;

return;
