home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / CPM / TURBOPAS / MAPSTATF.LBR / CLUSTER.PZS / CLUSTER.ÐAS
Text File  |  2000-06-30  |  8KB  |  278 lines
  1. (* Multivariate Analysis Package - Kmeans Clustering Module
  2.    Copyright 1986 Douglas L. Anderton.  This program may be freely
  3.    circulated so long as it is not sold for profit and any charge does
  4.    not exceed costs of reproduction *)
  5.  
  6. Program Cluster(Input,Output);
  7. Const
  8.    N=10; { number of vars in clustering }
  9.    C=25; { max-min number of clusters to consider, e.g. between 50 and
  10.            20 (or 100 and 70, etc.) clusters should be examined }
  11.    NBYC=250;
  12.    WORK=3500;
  13. Type
  14.    SUBS=1..N;
  15.    NCLS=1..C;
  16.    NSAV=1..NBYC;
  17.    RVEC = Array [SUBS] Of Real;
  18.    RVC2 = Array [NCLS] Of Real;
  19.    RVC3 = Array [NSAV] Of Real;
  20.    LVEC = Array [1..WORK] Of Real; { limited only by memory }
  21.    IVEC = Array [SUBS] Of Integer;
  22.    S8 = Array [SUBS] Of String[8];
  23. Var
  24.   dfile, ofile, tfile : Text;
  25.   sel : IVEC;
  26.   varn : S8;
  27.   vars, miss : RVEC;
  28.   i, j, k, l, nc, n1, n2, nv : Integer;
  29.   dv, ot, f : Byte;
  30.   x : LVEC;
  31.   s : RVC3;
  32.   e : RVC2;
  33.   d : Real;
  34.   df : String[1];
  35.   tfl : String[12];
  36.  
  37. {$I TRANSBUF.LIB}
  38.  
  39. Procedure openfiles(Var dfile, ofile:Text; Var ot:Byte);
  40. Var
  41.   dfl, ofl:String[12];
  42. Begin
  43.   ClrScr; Writeln(' *** CLUSTER: KMEANS CLUSTERING ***');
  44.   Writeln; Write('Name of the data file? ');
  45.   Readln(dfl); Assign(dfile,dfl); Reset(dfile);
  46.   Write('Name of the output file? ');
  47.   Readln(ofl); Assign(ofile,ofl); Rewrite(ofile);
  48.   ot:= 0;
  49.   If (ofl='CON:') Or (ofl='con:') Then ot:=1;
  50.   If (ofl='LST:') Or (ofl='lst:') Then ot:=2;
  51.   If (ot=2) Then
  52.     Begin
  53.     Writeln(ofile,'Multivariate Analysis Package (1.6) - ',
  54.       'Program: CLUSTER, Datafile: ',dfl); Writeln(ofile);
  55.     End;
  56.   End; (* Of openfiles *)
  57.  
  58. Procedure wait(ot:Byte);
  59. Begin
  60.   If ot=1 Then Begin
  61.     Write('- Press any key to continue -'); While Not KeyPressed Do; ClrScr;
  62.     End;
  63.   End; (* of wait *)
  64.  
  65. Procedure selectvar(Var sel:IVEC; Var varn:S8; Var miss:RVEC;
  66.                     Var nv:Integer; Var dv:Byte);
  67. Var
  68.    cfile:Text;
  69.    cfl:String[12];
  70.    i,j,f:Integer;
  71.    mis:Real;
  72.    van:String[8];
  73. Begin
  74.   Write('Name of the codebook file (or NONE)? '); Readln(cfl);
  75.   If (cfl<>'NONE') And (cfl<>'none') Then f:=1 Else f:=0;
  76.   If f=1 Then Begin Assign(cfile,cfl); Reset(cfile); End;
  77.   Writeln;
  78.   Write('How many variables in data file? '); Readln(nv);
  79.   Write('Number of variables to use in CLUSTER? '); Readln(dv);
  80.   Writeln('If multiple variables specified all should be comparable scale.');
  81.   For i := 1 To dv Do
  82.     Begin
  83.     Write('Column number for variable ',i,'? '); Readln(sel[i]);
  84.     Str(sel[i]:3,varn[i]); miss[i]:=-1E37;
  85.     End;
  86.    If f=1 Then Begin
  87.      For j:=1 to nv Do Begin
  88.        mis:=-1E37;
  89.        Readln(cfile,f,van,mis);
  90.        For i:=1 to dv Do
  91.          If f=sel[i] Then Begin
  92.            varn[i]:=van; miss[i]:=mis;
  93.            Writeln('Col: ',sel[i],'  Name: ',varn[i],' Missing: ',miss[i]:6);
  94.            End;     
  95.        End;
  96.      Close(cfile);
  97.      wait(1);
  98.      End;
  99.   End; (* Of selectvar *)
  100.  
  101. Procedure getcase(Var vars:RVEC; sel:IVEC; nv, dv:Integer; Var dfile:Text);
  102. Var
  103.   i, j:Integer;
  104.   x:Real;
  105. Begin
  106.   For i:=1 To nv Do
  107.     Begin
  108.     Read(dfile,x);
  109.     For j:=1 To dv Do If (sel[j]=i) Then vars[j]:=x;
  110.     End;
  111.   End; (* Of getcase *)
  112.  
  113. Procedure kmeans(Var x:LVEC; Var s:RVC3; Var e:RVC2; Var l, nc:Integer;
  114.                  Var dv:Byte; Var d:Real);
  115. Var
  116.   q: Array [NCLS] Of Integer;
  117.   i, it, j, k, r, v, w: Integer;
  118.   a, b, f, h, g, t: Real;
  119. Begin
  120.   For i:=1 To l Do
  121.     Begin
  122.     q[i]:=0; e[i]:=0.0;
  123.     For j:=1 To dv Do s[(i-1)*dv+j]:=0.0;
  124.     End;
  125.   For i:=1 To nc Do
  126.     Begin
  127.     k:=Trunc(x[(i-1)*(dv+1)+dv+1]+0.1); q[k]:=q[k]+1;
  128.     For j:=1 To dv Do s[(k-1)*dv+j]:=s[(k-1)*dv+j]+x[(i-1)*(dv+1)+j];
  129.     End;
  130.   For i:=1 To l Do
  131.     For j:=1 To dv Do s[(i-1)*dv+j]:=s[(i-1)*dv+j]/q[i];
  132.   d:=0.;
  133.   For i:=1 To nc Do
  134.     Begin
  135.     f:=0; k:=Trunc(x[(i-1)*(dv+1)+dv+1]+0.1);
  136.     For j:=1 To dv Do f:=f+Sqr(s[(k-1)*dv+j]-x[(i-1)*(dv+1)+j]);
  137.     e[k]:=e[k]+f;
  138.     d:=d+f;
  139.   End;
  140.   i:=0; it:=0;
  141.   While it<nc Do
  142.     Begin
  143.     i:=i+1; If i>nc Then i:=i-nc;
  144.     r:=Trunc(x[(i-1)*(dv+1)+dv+1]+0.1);
  145.     If q[r]>1 Then
  146.       Begin
  147.       h:=q[r]; h:=h/(h-1.0); f:=0;
  148.       For j:=1 To dv Do f:=f+Sqr(s[(r-1)*dv+j]-x[(i-1)*(dv+1)+j]);
  149.       a:=h*f; b:=1E30;
  150.       For j:=1 To l Do
  151.         Begin
  152.         If j<>r Then
  153.           Begin
  154.           h:=q[j]; h:=h/(h+1.0); f:=0;
  155.           For k:=1 To dv Do f:=f+Sqr(s[(j-1)*dv+k]-x[(i-1)*(dv+1)+k]);
  156.           f:=h*f;
  157.           If f<=b Then
  158.             Begin
  159.             b:=f; v:=j; w:=q[j];
  160.             End;
  161.           End;
  162.         End;
  163.       If b>=a Then
  164.         Begin
  165.         it:=it+1;
  166.         Write(' ',it);
  167.         End
  168.       Else
  169.         Begin
  170.         it:=0; e[r]:=e[r]-a; e[v]:=e[v]+b; d:=d-a+b;
  171.         h:=q[r]; g:=w; a:=1./(h-1.0); b:=1./(g+1.0);
  172.         For k:=1 To dv Do
  173.           Begin
  174.           s[(r-1)*dv+k]:=a*(h*s[(r-1)*dv+k]-x[(i-1)*(dv+1)+k]);
  175.           s[(v-1)*dv+k]:=b*(g*s[(v-1)*dv+k]+x[(i-1)*(dv+1)+k]);
  176.           End;
  177.         x[(i-1)*(dv+1)+dv+1]:=v; q[r]:=q[r]-1; q[v]:=q[v]+1;
  178.         End;
  179.       End;
  180.     End;
  181.   End; (* Of kmeans *)
  182.  
  183. Begin (* main *)
  184.   openfiles(dfile, ofile, ot);
  185.   selectvar(sel, varn, miss, nv, dv);
  186.   (* initialize *)
  187.   ClrScr; Writeln; n1:=0; n2:=C+1;
  188.   While (n2-n1)>C Do Begin
  189.     Write('Minimum number of Clusters to investigate? ');Read(n1);
  190.     Write(' Maximum number? ');Readln(n2);
  191.     If(n2-n1)>C Then Writeln('Too many, specify smaller range.');
  192.     End;
  193.   (* read in matrix *)
  194.   Writeln;
  195.   nc:=0;
  196.   While Not EOF(dfile) Do
  197.     Begin
  198.     getcase(vars, sel, nv, dv, dfile);
  199.     If Frac(nc/10)=0.0 Then Write('+');
  200.     If Not EOF(dfile) Then
  201.       Begin
  202.       f:=0;
  203.       For j:=1 To dv Do If vars[j]=miss[j] Then f:=1;
  204.       If f=0 Then
  205.         Begin
  206.         nc:=nc+1;
  207.         For j:=1 To dv Do
  208.           Begin
  209.           l:=(nc-1)*(dv+1)+j;
  210.           If l>WORK Then
  211.             Begin
  212.             Writeln('*** STORAGE EXCEEDED: Try Fewer (MAX-MIN) Clusters ***');
  213.             bdos(0);
  214.             End;
  215.           x[l]:=vars[j];
  216.           End;
  217.         End;
  218.       End;
  219.     End;
  220.   Reset(dfile);
  221.   (* compute And Output *)
  222.   For l:=n1 To n2 Do
  223.     Begin
  224.     k:=0;
  225.     For j:=1 To nc Do x[(j-1)*(dv+1)+dv+1]:=j-(j div l)*l+1;
  226.     ClrScr;
  227.     Writeln('Iterating on ',l:2,' clusters.');
  228.     kmeans(x,s,e,l,nc,dv,d);
  229.     ClrScr;
  230.     Writeln(ofile);
  231.     Writeln(ofile,'Total SSQ Distances Within ',l:2,' Clusters: ',d:10);
  232.     If ot<>1 Then
  233.       Writeln('Total SSQ Distances Within ',l:2,' Clusters: ',d:10);
  234.     wait(ot);
  235.     Writeln(ofile);
  236.     For i:=1 To l Do
  237.       Begin
  238.       Writeln(ofile,'Cluster: ',i:2,' Within SSQ: ',e[i]:10);
  239.       Writeln(ofile,' Centroids:');
  240.       For j:=1 To dv Do
  241.         Begin
  242.         Write(ofile,' ',sel[j]:2,': ',s[(i-1)*dv+j]:10);
  243.         If (j=5) And (dv>5) Then Writeln(ofile);
  244.         End;
  245.       Writeln(ofile);
  246.       If ((Frac(i/5)=0.0) Or (i=l)) Then wait(ot);
  247.       End;
  248.     Writeln('Would you like data and cluster membership for this level');
  249.     Write('written to a merged file (y/n)? '); Readln(df);
  250.     If (df='y') Or (df='Y') Then
  251.       Begin
  252.       Writeln;
  253.       Write('Name of new output file? '); Readln(tfl);
  254.       initwrite(tfl);
  255.       k:=0;
  256.       While Not EOF(dfile) Do
  257.         Begin
  258.         f:=0; For i:=1 To nv Do Read(dfile,s[i]);
  259.         If Frac(k/10)=0.0 Then Write('+');
  260.         If Not EOF(dfile) Then
  261.           Begin
  262.           For i:=1 To nv Do
  263.             For j:=1 To dv Do If (sel[j]=i) And (s[i]=miss[j]) Then f:=1;
  264.           s[nv+1]:=-1E37;
  265.           If f=0 Then
  266.             Begin k:=k+1; s[nv+1]:=x[(k-1)*(dv+1)+dv+1]; End;
  267.           For i:=1 To nv+1 Do Write(Usr,s[i]:9,' ');
  268.           Writeln(Usr);
  269.           End;
  270.         End;
  271.       Close(dfile); endwrite; Assign(dfile,tfl); Reset(dfile); nv:=nv+1;
  272.       End;
  273.     ClrScr;
  274.     End;
  275.   Close(dfile); Close(ofile);
  276.   Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
  277. End.
  278. wait(ot)