proc format;
invalue animal '+'=1 '-'=0;            
DATA GEL;
input  id $ 1-2  popul $ 3-4 
       @5 (band1-band10) (animal.);      
if mean(of band1-band10)= . then delete; 

cards;
A Y +++++++++-
B Y +++++++++-
C X ----++---+
D Y ++-+++--+-
E X ++--++---+
F Y +++++++++-
G Y ++-+++--+-
H Y +++++++++-
I Y ++++++-++-
J X ----++---+
K Y +++++++++-
L X --+++----+
;
run;

proc sort data=gel out=gel; by popul id;  
proc print data=gel noobs;
                                

proc summary;                   
class popul;
var band1;
output out=nn n=nn;

DATA GEL;                      
SET GEL;                       
*if xxx then delete;
DROP POPUL ID;

PROC IML;
USE GEL;
READ ALL INTO A;

use nn;
read all into nn;
na=nn[2,3];
nb=nn[3,3];
nab=na+nb;

print nn na nb;
   C=(A*A`);                   
   PRINT C[format=2.0];

   N = NCOL(C);
   D = diag(C);
   O = C-D;
   rowi=J(1, N, 1);

   SUMC= rowi * O * rowi`;
   SUMD= rowi * D * rowi`;

     BS=J(N,N,0);
     do i= 1 to N  by 1;
        do j= (i+1) to N  by 1;
          BS[i,j] =(2*C[i,j])/ (C[i, i] + C[j, j]);
        end;
    end;

   BS_SUM   = (SUMC / (SUMD*(N-1)));       
   BS_Lynch = 2*(rowi * BS * rowi`) /(N*(N-1));

   PRINT BS[format=5.3];   
   PRINT BS_SUM BS_Lynch;   
   

                            
   rowa=J(1,na,1)||J(1,nb,0);
   rowb=J(1,na,0)||J(1,nb,1);                 
   BSx= 2*(rowa * BS * rowa`) /(na*(na-1));   
   BSy= 2*(rowb * BS * rowb`) /(nb*(nb-1));  
   BSxy=  (rowa * BS * rowb`) /(na*nb);        

   B_Lynch= 1 + BSxy - 0.5*(BSx+ BSy);         
   Dxy = -log (BSxy /SQRT (BSx * BSy) );       
   print BSx BSy BSxy na nb B_Lynch Dxy;

   D = (A`*A);                  
   PRINT D[format=2.0];         
   M = NCOL(D);                 
   B = J(M, N, 1);
   B = (B - A`);
   E = (B*B`);                  
   F = (A`*B`);

   E=E+J(M, M, 1);             
                                

   R = J(M, M, 0);             
     do i= 1 to M  by 1;
        do j= 1 to M  by 1;
         B=(D[i, j]*E[i, j] - F[i, j]*F[j, i]);
         K=SQRT( (D[i, j]+F[i,j])*(F[j,i]+E[i,j])
               *(D[i, j]+F[j,i])*(F[i,j]+E[i,j]));
                  if K<>0 then  R[i,j] = B/K;
              if((SQRT( (1.0 - R[i,j]*R[i,j])/(N-2))) <> 0) 
                then T=R[i,j]/(SQRT( (1.0 - R[i,j]*R[i,j])/(N-2)));
           if(R[i,j]=1.0)                 then T=999;
 *        if(PROBT(abs(T), N-2) < 0.975) then R[i,j]=0.000; 
        end;                             
      end;

  PRINT R[format=5.3];
  PRINT D[format=4.0] E[format=4.0] F[format=4.0];

  R= J(M, M, 1) - abs(R);       
                                0=>1, 0.1=>0.9, -0.9=>0.1, 1=>0 
  create corr from R;           
  append from R;

  H=A`;                         
  create BAND from H;          
  append from H;

PROC CLUSTER DATA=corr(type=distance) METHOD=SINGLE;
PROC TREE HORIZONTAL SPACES=2; 

PROC CLUSTER DATA=BAND METHOD=AVERAGE;
PROC TREE HORIZONTAL SPACES=2;
RUN;