makeevenCIJ.m 1.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475
  1. function [CIJ] = makeevenCIJ(N,K,sz_cl)
  2. %MAKEEVENCIJ Synthetic modular small-world network
  3. %
  4. % CIJ = makeevenCIJ(N,K,sz_cl);
  5. %
  6. % This function generates a random, directed network with a specified
  7. % number of fully connected modules linked together by evenly distributed
  8. % remaining random connections.
  9. %
  10. % Inputs: N, number of vertices (must be power of 2)
  11. % K, number of edges
  12. % sz_cl, size of clusters (power of 2)
  13. %
  14. % Outputs: CIJ, connection matrix
  15. %
  16. % Notes: N must be a power of 2.
  17. % A warning is generated if all modules contain more edges than K.
  18. % Cluster size is 2^sz_cl;
  19. %
  20. %
  21. % Olaf Sporns, Indiana University, 2005/2007
  22. % compute number of hierarchical levels and adjust cluster size
  23. mx_lvl = floor(log2(N));
  24. sz_cl = sz_cl-1;
  25. % make a stupid little template
  26. t = ones(2).*2;
  27. % check N against number of levels
  28. Nlvl = 2^mx_lvl;
  29. if (Nlvl~=N)
  30. disp('Warning: N must be a power of 2');
  31. end;
  32. N = Nlvl;
  33. % create hierarchical template
  34. for lvl=1:mx_lvl-1
  35. CIJ = ones(2^(lvl+1),2^(lvl+1));
  36. group1 = 1:size(CIJ,1)/2;
  37. group2 = size(CIJ,1)/2+1:size(CIJ,1);
  38. CIJ(group1,group1) = t;
  39. CIJ(group2,group2) = t;
  40. CIJ = CIJ+ones(size(CIJ,1),size(CIJ,1));
  41. t = CIJ;
  42. end;
  43. s = size(CIJ,1);
  44. CIJ = CIJ-ones(s,s)-mx_lvl.*eye(s);
  45. % assign connection probabilities
  46. %CIJp = mx_lvl-CIJ-sz_cl;
  47. %CIJp = (CIJp>0).*CIJp;
  48. CIJp = (CIJ>=(mx_lvl-sz_cl));
  49. % determine number of remaining (non-cluster) connections and their
  50. % possible positions
  51. %CIJc = (CIJp==0);
  52. CIJc = (CIJp==1);
  53. remK = K-nnz(CIJc);
  54. if (remK<0)
  55. disp('Warning: K is too small, output matrix contains clusters only');
  56. end;
  57. [a,b] = find(~(CIJc+eye(N)));
  58. % assign 'remK' randomly distributed connections
  59. rp = randperm(length(a));
  60. a = a(rp(1:remK));
  61. b = b(rp(1:remK));
  62. for i=1:remK
  63. CIJc(a(i),b(i)) = 1;
  64. end;
  65. % prepare for output
  66. CIJ = CIJc;