%**************************************************;
%* Olivier Godechot, V1.5 le 19.04.2008           *;
%**************************************************;
%* Macro calculating Burt s structural constraint.*;
%**************************************************;


%Macro burtcstr(dataset=,orig=,dest=,value=,form=Yes);

%window Welcome
#1 @1 "This macro calculates the burt constraint for a network"
#3 @1 "Data must be organized the following way : "
#5 @1 "Origin_id 	Destiny_id 	Value"
#6 @1 "a1572 		a1759 			2"
#7 @1 "a1088 		a2716 			1"
#8 @1 "a1863 		a2716 			1"
#10 @1 "Reading : An arc of value 2 is going from vertex a1572 to vertex a1759"
#12 @1 "What is the name of your SAS dataset ? (data=)" @75 dataset 30 attr=rev_video auto=yes
#14 @1 "What is the variable name that identifies the origin of an arc ? (orig=)" @75 orig 30 attr=rev_video auto=yes
#16 @1 "What is the variable name that identifies the destiny of an arc ? (dest=)" @75 dest 30 attr=rev_video auto=yes
#18 @1 "What is the variable name for the value of an arc ? (value=)" @75 value 30 attr=rev_video auto=yes

#20 @1 "Once completed, press Enter "

#22 @1 %nrstr("[Complete syntax: %burtcstr(dataset=,orig=,dest=,value=,form=); ]")
#23 @1 "* Required. Id variables must not contain the character  # "
#30 @50 "Credits : Olivier GODECHOT. http://olivier.godechot.free.fr/ "
;

%if &form=Yes %then %display Welcome;

%let dat0=_%scan(&dataset,1)_C0;
%let dat1=_%scan(&dataset,1)_C1;
%let dat2=_%scan(&dataset,1)_Ci;
%let dat3=_%scan(&dataset,1)_C3;
%let dat4=_%scan(&dataset,1)_C4;
%let dat5=_%scan(&dataset,1)_cij;

%* let orig_=%str(&orig)_;
%let dest_=%scan(&dest,1);

proc iml;
use &dataset;
  read all var{&orig} into orig;
  read all var{&dest} into dest;

%if &value= %then %do;
	value=J(nrow(orig),1,1);
%end;

%else  %do;
   read all var{&value} into value;
%end;

%* adj fonction from J. Moody;
start adj(snd,rcv,value);
   nomset=unique(snd,rcv);
   if type(nomset)='C' then do;
      nomset=setdif(nomset,'.');
   end;
   else do;
     nomset=setdif(nomset,.);
   end;
   adjmat=j(ncol(nomset),ncol(nomset),0);
   do i=1 to nrow(snd);
    sendloc=loc(nomset=snd[i]);
    if type(sendloc)='N' then do;
      rcvset=unique(rcv[i,]);
      if type(rcvset)='C' then do;
         rcvset=setdif(rcvset,'.');
      end;
      else do;
       rcvset=setdif(rcvset,.);
      end;
       if type(rcvset)^='U' then do;
          do j=1 to ncol(rcvset);
           jloc=loc(nomset=rcvset[j]);
           adjmat[sendloc,jloc]=adjmat[sendloc,jloc]+value[i];
          end;
       end;
     end;
    end;

 nomset=nomset`;
 if type(nomset)='N' then do;
    adjmat=nomset||adjmat;
 end;
 else do;
   print 'Character values can not be appended to Adjacency Matrix.';
   print 'The nodes are thus labeled from 1 to g, in the following order:';
   idx=1:nrow(adjmat);
   idx=idx`;
   print nomset idx;
   adjmat=idx||adjmat;
 end;
 return(adjmat);
finish;
%* reach fonction from J. Moody;
start reach (inmat);
 r=inmat;
 rt=r;
 t=2;
 do until (sdif=0);
  rt=rt*inmat;  /* multiply tmat by inmat=taking it to the next power */
  rt0=rt-diag(rt); /* give matrix with diag=0 */
  rt01=rt0>0;   /* make nonzero elements=1 */
  mark=rt01>r;
  tmark=t#mark; /* makes a replacement matrix of vavlue t */
  k=tmark+r;
  sk=sum(k);
  sr=sum(r);
  sdif=sum(sk-sr);  /* when this is zero, no new paths can be made */
  t=t+1;
  r=tmark+r;
  free drt rt0 rt01 mark k sk sr;
 end;
 return(r);
finish;

%* My fonction;
start matid(id);
n=nrow(id);
do i=1 to nrow(id);
	do j=1 to nrow(id);
	couple=id[i] ||	id[j];
		matid=matid || T(couple);
	end;
end;
return(T(matid));
finish;


  adjmat=adj(orig,dest,value);
  adjmat=adjmat[,2:ncol(adjmat)];

	adjmat=adjmat+T(adjmat);



* Size of the matrice;
n=nrow(adjmat);
matb=adjmat#(J(n,n,1)-I(n));

* Investment rate in one given contact;
invest=(matb+T(matb))#T(shape(((matb*shape(1,n,1)+T(matb)*shape(1,n,1))##-1),n,n));

symdich=invest>J(n,n,0);

* Importance of the constraint from j towards vers i ;
Cij=((invest+invest*invest)##2)#symdich;

id=T(unique(orig,dest));

* Sum of constraints;
C= cij*J(n,1,1);
print C;

*Creation of a sas dataset;
cijbis= T(shape(cij,1));
matid=matid(id);



create &dat0 from id[colname={ &orig }];
append from id;

create &dat1 from C[colname={constraint}];
append from C;

create &dat3 from matid[colname={ &orig &dest_}];
append from matid;

create &dat4 from cijbis[colname={ cij}];
append from cijbis;

quit;

data &dat2; merge &dat0 &dat1;
run;

data &dat5; merge &dat3 &dat4;
run;


proc means data=&dat2 N NMISS MEAN STD MIN Q1 MEDIAN Q3 MAX ; var constraint;
Title "Statistics on structural constraints ";
run;

DATA _NULL_;
file print;
PUT "Structural contraints are stored in the &dat2 dataset.";
PUT "Contribution of j on i� s structural contraint are stored in the &dat1 dataset.";
PUT "The macro program was ";
put %nrstr("%burtcstr");
PUT "(dataset=&dataset,orig=&orig,dest=&dest,value=&value,form=&form);";
run;

proc datasets; delete &dat0 &dat1 &dat3 &dat4; quit;

%MEND;