(* Multivariate Analysis Package - Principle Component Factoring Module
   Copyright 1985 Douglas L. Anderton.  This program may be freely
   circulated so long as it is not sold for profit and any charge does
   not exceed costs of reproduction *)

Program Factor(Input,Output);
Const
  N=30;
Type
  SUBS=1..N;
  RVEC = Array [SUBS] Of Real;
  NBYN = Array [SUBS] Of RVEC;
  IVEC = Array [SUBS] Of Integer;
  S8 = Array [SUBS] Of String[8];
Var
  dfile, ofile : Text;
  sel : IVEC;
  e, u, v, w, x, y: RVEC;
  varn : S8;
  cor, a : NBYN;
  nf, nv, dv, i, j, k, rn, ot: Integer;
  tol, temp, tmp2, nc : Real;
  error : boolean;

Procedure openfiles(Var dfile, ofile:Text; Var nv, dv, ot:Integer);
Var
  dfl, ofl:String[12];
Begin
  ClrScr; Writeln(' *** FACTOR: PRINCIPLE COMPONENT FACTORING ***');
  Writeln; Write('Name of the CORREL file? ');
  Readln(dfl); Assign(dfile,dfl); Reset(dfile);
  Write('Name of the output file? ');
  Readln(ofl); Assign(ofile,ofl); Rewrite(ofile);
  ot:= 0;
  If (ofl='CON:') Or (ofl='con:') Then ot:=1;
  If (ofl='LST:') Or (ofl='lst:') Then ot:=2;
  If (ot=2) Then
    Begin
    Writeln(ofile,'Multivariate Analysis Package (1.6) - ',
      'Program: FACTOR, Datafile: ',dfl); Writeln(ofile);
    End;
  Write('How many variables in CORREL matrix? '); Readln(nv);
  Write('Number of variables to use in FACTOR? '); Readln(dv);
  End; (* Of openfiles *)

Procedure selectvar2(Var sel:IVEC; dv:Integer);
Var
  i:Integer;
Begin
  Writeln;
  For i:=1 To dv Do Begin
    Write('Column number for variable ',i,'? '); Readln(sel[i]); End;
  End; (* Of selectvar *)

Procedure wait(ot:Integer);
Begin
  If ot=1 Then Begin
    Write('- Press any key to continue -'); While Not KeyPressed Do; ClrScr;
    End;
  End; (* of wait *)

{$I GETCOR.LIB}
{$I EIGEN.LIB}
{$I MATINV.LIB}

Procedure rotate(Var a:NBYN; nf, nv, meth, ot:Integer);
Var
  c:RVEC;
  i,j,k,f,nr:Integer;
  ep,em,a1,b1,c1,d1,u,v,qn,qd,cs,sn,cs1,sn1,sp,cp:Real;
Begin
  Writeln(ofile);
  If(meth=0) Then Writeln(ofile,'Varimax Rotation of ',nf:2,' Factors')
  Else Writeln(ofile,'Quartimax Rotation of ',nf:2,' Factors');
  ep:=0.00116;
  For j:=1 To nv Do
    Begin
    c[j]:=0.0;
    For k:=1 To nf Do c[j]:=c[j]+Sqr(a[j,k]);
    c[j]:=Sqrt(c[j]);
    For k:=1 To nf Do a[j,k]:=a[j,k]/c[j];
  End;
  While((nr-((nf*(nf-1))/2))<>0) Do
    Begin
    nr:=0;
    For i:=1 To (nf-1) Do
      Begin
      For j:=(i+1) To nf Do
        Begin
        a1:=0.0; b1:=0.0; c1:=0.0; d1:=0.0;
        For k:=1 to nv Do Begin
          u:=Sqr(a[k,i])-Sqr(a[k,j]); v:=a[k,i]*a[k,j]*2.0;
          a1:=a1+u; b1:=b1+v; c1:=c1+Sqr(u)-Sqr(v); d1:=d1+(u*v*2.0); End;
        If(meth=0) Then
          Begin qn:=d1-((2*a1*b1)/nv); qd:=c1-((Sqr(a1)-Sqr(b1))/nv); End;
        If(meth=1) Then Begin qn:=d1; qd:=c1; End;
        f:=0;
        If((Abs(qn)+Abs(qd))=0) Then Begin nr:=nr+1; f:=1; End;
        If((Abs(qn)-Abs(qd))=0) and (f=0) Then
          Begin cs:=0.70710678; sn:=cs; f:=2; End;
        If((Abs(qn)-Abs(qd))>0) and (f=0) Then
          Begin
          em:=Abs(qd/qn);
          If(em<ep) Then Begin cs:=0; sn:=1; End
          Else Begin sn:=1/Sqrt(1+Sqr(em)); cs:=sn*em; End;
          f:=2;
          End;
        If(f=0) Then
          Begin
          em:=Abs(qn/qd);
          If(ep<=em) Then Begin
            cs:=Cos(ArcTan(em)); sn:=Sin(ArcTan(em)); f:=2; End
          Else If(qd>=0) Then Begin nr:=nr+1; f:=1; End
            Else Begin sp:=0.70710678; cp:=sp; f:=3; End
          End;
        If(f=2) Then
          Begin
          em:=Sqrt((1+cs)*0.5); cs1:=Sqrt((1+em)*0.5); sn1:=sn/(4*cs1*em);
          If(qd<0) Then
            Begin cp:=0.70710678*(cs1+sn1); sp:=0.70710678*(cs1-sn1); End
          Else Begin cp:=cs1; sp:=sn1; End;
          If(qn<0) Then sp:=-sp;
          End;
        If((f=2) or (f=3)) Then
          For k:=1 To nv Do
            Begin em:=(a[k,i]*cp)+(a[k,j]*sp);
            a[k,j]:=(a[k,j]*cp)-(a[k,i]*sp); a[k,i]:=em; End;
        End;
      End;
    End;
  For i:=1 To nv Do For j:=1 To nf Do a[i,j]:=a[i,j]*c[i];
  For j:=1 To nf Do
    Begin c[j]:=0.0; For k:=1 To nv Do c[j]:=c[j]+Sqr(a[k,j]); End;
  ClrScr;
  For i:=1 To nf Do Begin
    Writeln(ofile); Writeln(ofile,'Rotated Factor Pattern For Factor: ',i:2);
    For j:=1 To dv Do
      Begin Writeln(ofile,'   ',varn[j]:8,': ',a[j,i]:10:8); End;
    Writeln(ofile,'Factor Accounts For ',(c[i]*100/nv):8:5,'% of Variance');
    wait(ot);
    End;
  End; (* of rotate *)

Begin (* main *)
  openfiles(dfile, ofile, nv, dv, ot);
  Write('Number of Factors to Extract (<=number of variables)? ');
  Readln(nf); If (nf=dv) Then rn:=0 Else rn:=1;
  tol:=1E-8;
  Write('Tolerance for eigenvalue convergence <return for 1E-8>? ');
  Readln(tol);
  selectvar2(sel, dv);
  (* read correlation matrix matrix *)
  getcor(cor, x, y, varn, nc, sel, nv, dv, dfile);
  (* full rank analysis *)
  If(rn=0) Then
    Begin
    a:=cor;
    (* invert *)
    Writeln('full rank analysis inverting matrix...');
    matinv(a,dv,temp,error);
    ClrScr; Writeln(ofile,'Inverse:  (determinant=',temp,')');
    If(error) Then Writeln(ofile,'Error reported in solving for Inverse');
    temp:=-((nc-1.0)-((1.0/6.0)*(2.0*nc+5.0)))*Ln(temp);
    Write(ofile,'Chi-Square for Sphericity = ',temp:10:3);
    Writeln(ofile,'  degrees of freedom = ',((Sqr(dv)-dv)/2):6:0);
    wait(ot);
    For j:=1 To dv Do w[j]:=1.0-(1.0/a[j,j]);
    End;
  a:=cor;
  eigen(cor, a, e, tol, dv, nf, ot);
  wait(ot);
  temp:=0.0; tol:=dv;
  For j:=1 To nf Do Begin
    u[j]:=(e[j]/tol)*100.0; temp:=temp+u[j]; v[j]:=temp;
    tmp2:=Sqrt(e[j]);
    For k:=1 To dv Do
      Begin cor[k,j]:=a[k,j]*tmp2; a[k,j]:=a[k,j]/tmp2; End;
    End;
  ClrScr; Writeln(ofile);
  Writeln(ofile,'Factor        Eigen          %Variance     %Cumulative');
  Writeln(ofile);
  For i:=1 To nf Do Writeln(ofile,'   ',i:2,'         ',e[i]:8:5,
   '       ',u[i]:7:3,'       ',v[i]:7:3);
  Writeln(ofile);
  wait(ot);
  For i:=1 To nf Do Begin
    Writeln(ofile); Writeln(ofile,'Factor Pattern For Factor: ',i:2);
    For j:=1 To dv Do
      Begin Writeln(ofile,'   ',varn[j]:8,': ',cor[j,i]:10:8); End;
    wait(ot);
    End;
  Writeln(ofile); Write(ofile,'Variable      Communality');
  If(rn=0) Then Writeln(ofile,'     R-Square') Else Writeln(ofile);
   For j:=1 To dv Do
    Begin temp:=0.0;
    For k:=1 To nf Do temp:=temp+Sqr(cor[j,k]);
{DOUG    If temp>1.0 then temp:=1.0; }
    Write(ofile,'   ',varn[j]:8,'    ',temp:10:7);
    If(rn=0) Then Writeln(ofile,'     ',w[j]:10:7) Else Writeln(ofile);
    End;
  wait(ot);
  For i:=1 To nf Do Begin
    Writeln(ofile); Writeln(ofile,'Score Coefficients For Factor: ',i:2);
     For j:=1 To dv Do
      Begin Writeln(ofile,'   ',varn[j]:8,': ',a[j,i]:10:8); End;
    wait(ot);
    End;
  If(nf>1) Then
    Begin
    Write('What type of Factor Rotation (0=none,1=varimax,2=quartimax)?');
    Readln(i); i:=i-1; If(i>=0) Then rotate(cor, nf, dv, i, ot);
    End;
   Close(ofile); Close(dfile);
   Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
  End.