function a=ws(n,r,p)
  for i=1:n-r for j=i+1:i+r a(i,j)=1; a(j,i)=1; end; end
  for i=n-r+1:n
    jind=[i+1:n, 1:i+r-n]
    for j=jind
      a(i,j)=1; a(j,i)=1; 
    end
  end
  v=rand(1,2*n*r); k=1
  for i=1:n-r
   j1ind=[1:i-1, i+1:n] 
   for j=i+1:i+r
     if v(k)<p then 
       m=0; ind=[]
       for j1=j1ind if a(i,j1)==0 ind(m+1)=j1; m=m+1; end; end 
       mr=int(v(k+1)*m)+1; jn=ind(mr); a(i,j)=0; a(j,i)=0; a(i,jn)=1; a(jn,i)=1
     end
     k=k+2;
   end
  end
  for i=n-r+1:n
    jind=[i+1:n, 1:i+r-n]; j1ind=[1:i-1, i+1:n]
   for j=jind
     if v(k)<p then 
       m=0; ind=[]
       for j1=j1ind if a(i,j1)==0 ind(m+1)=j1; m=m+1; end; end 
       mr=int(v(k+1)*m)+1; jn=ind(mr); a(i,j)=0; a(j,i)=0; a(i,jn)=1; a(jn,i)=1
     end
     k=k+2;
   end
  end
endfunction

function netplot(a)
  n=size(a,"r"); th=2*%pi/n*(1:n); z=exp(%i*th); xbasc(); drawlater()
  for i=1:n-1
    for j=i+1:n
      if a(i,j)==1 plot([real(z(i)), real(z(j))], [imag(z(i)), imag(z(j))]); end
    end
  end
  drawnow()
endfunction

function d=deg(a)
  n=size(a,"r")
  for i=1:n d(i)=sum(a(i,:)); end
endfunction

function c=cluster(a)
  n=size(a,"r"); v=[]
  for i=1:n
   m=1  
   for j=1:n 
    if a(i,j)==1 v(m)=j; m=m+1; end
   end
   b=a(v,v); k=m-1
   if k>1 then c(i)=sum(b)/k/(k-1); else c(i)=0; end
   v=[]
  end
endfunction

function [l, nk]=dist(a,lmax)
  n=size(a,"r"); l=zeros(1,n*(n-1)/2); ak=a
  k=1; flag=%t; nk=[]
  while flag&(k<=lmax)
    m=0; s=0;
    for i=1:n
      for j=i+1:n
        s=s+1;
        if (ak(i,j)>0)&(l(s)==0) then l(s)=k; m=m+1; end
      end
    end
    if m>0 then nk(k)= m; ak=ak*a; k=k+1; else flag=%f; end
  end
endfunction

stacksize(10000000)
ini=["00100";"2";"0.1000"];    // valori assegnati per default
[ok,n,r,p]=getvalue("Dare i valori dei parametri",..
  ["n";"r";"p"],list("col",1,"col",1,"col",1),ini);

msprintf(' n = %f,  r= %f\n', n, r)
a=ws(n,r,p);
mean_deg=sum(a)/n;
mprintf(' grado medio = %f\n', mean_deg)
ncl=0; d=deg(a);
for i=1:n if d(i)>=2 ncl=ncl+1; end; end
cluster_coeff=sum(cluster(a))/ncl;
mprintf(' coeff. di clustering = %f\n', cluster_coeff)
[l nk]=dist(a,20);
mean_dist=sum(l)/sum(nk);
mprintf(' distanza media = %f\n', mean_dist)
conn_coeff=2*sum(nk)/n/(n-1);
mprintf(' perc. di coppie connesse = %f\n', conn_coeff)
netplot(a)
