(* Multivariate Analysis Package - Correlation and Covariance 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 Correl(Input,Output);
Const
   N=50;
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;
  cvj, muj, nc, miss, vars, mean, stdev, sm : RVEC;
  varn : S8;
  mu, cv : NBYN;
  nv, i, j, k, wgt : Integer;
  ot, dv, wt : SUBS;
  vj, vw, vr, wg2, wht, corr : Real;
  wflg : Boolean;

Procedure openfiles(Var dfile, ofile:Text; Var ot:SUBS);
Var
  dfl, ofl:String[12];
Begin
  ClrScr; Writeln(' *** CORREL: CORRELATION AND COVARIANCE MATRIX ***');
  Writeln; Write('Name of the data 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: CORREL, Datafile: ',dfl); Writeln(ofile);
    End;
  End; (* Of openfiles *)

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 *)

Procedure selcvar(Var sel:IVEC;Var varn:S8;Var miss:RVEC;
                  Var nv:Integer; Var dv, wt:SUBS);
Var
  cfile:Text;
  cfl:String[12];
  i,j,f:Integer;
  mis:Real;
  van:String[8];
Begin
  Write('Name of the codebook file (or NONE)? '); Readln(cfl);
  If (cfl<>'NONE') And (cfl<>'none') Then f:=1 Else f:=0;
  If f=1 Then Begin Assign(cfile,cfl); Reset(cfile); End;
  Writeln;
  Write('How many variables in data file? '); Readln(nv);
  Write('Number of variables to use in CORREL? '); Readln(dv);
  For i := 1 To dv Do
    Begin
    Write('Column number for variable ',i,'? '); Readln(sel[i]);
    Str(sel[i]:3,varn[i]); miss[i]:=-1E37;
    End;
  Write('Of these Column numbers which is weight (0=none)? '); Readln(wt);
  If f=1 Then Begin
    For j:=1 to nv Do Begin
      mis:=-1E37;
      Readln(cfile,f,van,mis);
      For i:=1 to dv Do
        If f=sel[i] Then Begin
          varn[i]:=van; miss[i]:=mis;
          Writeln('Col: ',sel[i],'  Name: ',varn[i],' Missing: ',miss[i]:6);
          End;
      End;
    Close(cfile);
    wait(1);
    End;
  End; (* Of selcvar *)

Procedure getcase(Var vars:RVEC;Var sel:IVEC; nv:Integer; dv:SUBS;
                  Var dfile:Text);
Var
  i, j:SUBS;
  x:Real;
Begin
  For i := 1 To nv Do
    Begin
    Read(dfile,x); For j:=1 to dv Do If (sel[j]=i) Then vars[j]:=x;
    End;
  End; (* Of getcase *)

Begin (* main *)
  openfiles(dfile, ofile, ot);
  selcvar(sel, varn, miss, nv, dv, wt);
  (* initialize *)
  wgt:=0; wflg:=false;
  FillChar(nc,6*N,0); FillChar(stdev,6*N,0); FillChar(sm,6*N,0);
  FillChar(cv,6*N*N,0); FillChar(mu,6*N*N,0);
  For i := 1 To dv Do
    If wt=sel[i] Then Begin wgt:=i; wflg:=true; End;
  (* accumulate sscp *)
  Writeln;
  k:=0;

  While Not EOF(dfile) Do
    Begin
    k:=k+1;
    getcase(vars, sel, nv, dv, dfile);
    If Frac(k/10)=0.0 Then Write('+');
    { core computational loop }
    vw:=vars[wgt]; { min subscript refs }
    If Not EOF(dfile) Then
      If Not (wflg And (vw=miss[wgt])) Then
      For j:=1 To dv Do
        If vars[j]<>miss[j] Then
          Begin
          vj:=vars[j]; { min subscript refs }
          If wflg And (j<>wgt) Then wht:=vw Else wht:=1.0;
          { spicer algorithm for descriptive stats }
          nc[j]:=nc[j]+wht;
          If wflg Then Begin { avoid mult and div by 1 if no weight }
            vr:=(vj-mean[j])*wht/nc[j];
            stdev[j]:=stdev[j]+(nc[j]*(nc[j]-wht)/wht)*Sqr(vr);
          End Else Begin
            vr:=(vj-mean[j])/nc[j];
            stdev[j]:=stdev[j]+(nc[j]*(nc[j]-1.0))*Sqr(vr);
          End;
          mean[j]:=mean[j]+vr;
          { loop for covariance and pairwise deletion means stored
          in same matrix as triangulars plus diagonal sm }
          cvj:=cv[j]; muj:=mu[j]; {slicing to avoid subs nc*dv*dv*6 times}
          For i:=j To dv Do
            If vars[i]<>miss[i] Then
              Begin
              If wflg And (i<>wgt) Then wg2:=wht*vw Else wg2:=wht;
              If i=j Then Begin sm[i]:=sm[i]+wg2; vr:=wg2/sm[i]; End
              Else {i>j}
                Begin
                cv[i,j]:=cv[i,j]+wg2; vr:=wg2/cv[i,j];
                mu[i,j]:=mu[i,j]+(vars[i]-mu[i,j])*vr;
                End;
              muj[i]:=muj[i]+(vj-muj[i])*vr;
              cvj[i]:=cvj[i]+(vj*vars[i]-cvj[i])*vr;
              End;
          cv[j]:=cvj; mu[j]:=muj;
          End;
    End;

  (* compute And Output *)
  ClrScr;
  Writeln(ofile,'Means and Standard Deviations:');
  Writeln(ofile);
  For j:=1 To dv Do
    Begin
    If nc[j] > 1.0 Then
      Begin
      stdev[j]:=Sqrt(stdev[j]/(nc[j]-1.0));
      Writeln(ofile,varn[j]:8,'  Mean:',mean[j]:13:5,
       '  Standard Deviation:',stdev[j]:13:5,'  Cases:',nc[j]:9:0);
      End;
    End;
  Writeln(ofile);
  wait(ot);
  Writeln(ofile,'Covariance and Correlations (pairwise deletion):');
  Writeln(ofile);
  For j:= 1 To dv Do
    If sm[j] > 1 Then
      cv[j,j]:=cv[j,j]-Sqr(mu[j,j]);
  For j:= 1 To dv Do
    Begin
    cvj:=cv[j]; muj:=mu[j]; {slicing}
    For i:=j+1 To dv Do
      If cv[i,j]>0 Then cvj[i]:=cvj[i]-(muj[i]*mu[i,j]);
    For i:=j+1 To dv Do
      Begin
      If cv[i,j]>0.0 Then
        Begin
        corr:=cvj[i]/Sqrt(cv[i,i]*cvj[j]);
        Writeln(ofile,varn[j]:8,' with ',varn[i]:8,'   Covariance: ',
          cvj[i]:13,'   Correlation:',corr:10:6);
        End;
      End;
    End;
  Close(dfile); Close(ofile);
  Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
End.