nsw_clusterSW_sel_link2

   SynGrid
   Copyright (c) 2008, 2017, Electric Power and Energy Systems (EPES) Research Lab
   by Zhifang Wang, Virginia Commonwealth University

   This file is part of SynGrid.
   Covered by the 3-clause BSD License (see LICENSE file for details).

0001 function [Ma] = nsw_clusterSW_sel_link2(Mda,degModel,mdeg)
0002 
0003 %   SynGrid
0004 %   Copyright (c) 2008, 2017, Electric Power and Energy Systems (EPES) Research Lab
0005 %   by Zhifang Wang, Virginia Commonwealth University
0006 %
0007 %   This file is part of SynGrid.
0008 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0009 
0010 % output "Ma"  is upper triangle matrix.
0011 N = length(Mda(:,1));
0012 Ma  = zeros(N); % --> adjacency matrix according to Nb geometric distribution of branch numbers
0013 if(degModel ==0)
0014     Nb = (geornd(1/(mdeg-1),[1 N])+1);%--> expected number of links per node, geometric distribution,
0015 elseif(degModel == 1)
0016     Nb = shifted_geornd(mdeg,N);
0017 elseif(degModel ==2)
0018 %     Nb = (poissrnd(mdeg-1,[1 N])+1);%--> expected number of links per node, Poisson distribution
0019     Nb = (sg_poissrnd(mdeg-1,1,N)+1);%--> expected number of links per node, Poisson distribution
0020 elseif(degModel ==3)
0021 else
0022     disp('wrong setting of p_brh, must be 1,2 or 3')
0023     return;
0024 end
0025 
0026 if(degModel==3)
0027     Ma = triu(Mda,1); % level-1 upper triangle sub-matrix,contains all connections in Mda
0028 else
0029     for ni = 1:N-1
0030         Ny = Nb(ni); % --> number of needed branches for node-ni
0031         if(ni>1) Ny = Ny - sum(Ma(1:ni-1,ni )); end %--> mimus previous branches with nodes (1,...,ni-1)
0032         if(Ny<=0)
0033             Ma(ni,ni+1:N)=zeros(1,N-ni);
0034             continue;
0035         end
0036 
0037         pos = ni+ find(Mda(ni,ni+1:N)>0); % available branches with nodes (ni+1,...,N)
0038         Na = length(pos);
0039         if(Na == 0 || Na <= Ny)
0040             %disp(['Na Ny',num2str(Na), '   ', num2str(Ny)]);
0041             Ma(ni,ni+1:N) = Mda(ni,ni+1:N);
0042             continue;
0043         end
0044         if(Na > Ny)
0045             %pos_Ny = randsrc(1,Ny,[pos, 1/Na*ones(1,Na)]);
0046             pos_Ny = randi([1 Na],1,Ny);
0047             Ma(ni,pos(pos_Ny)) = 1;
0048         end
0049     end % --> for ni = 1:N-1
0050 end % --> if(degModel==3)
0051 
0052 
0053 %--------------------------------------------------------------------------
0054 function yout = shifted_geornd(mdeg, N)
0055 % shifted Geometric distribution
0056 % p1 = 0.1, for degree = 1 nodes;  <-- group1
0057 % p2 = 1-p1 = 0.9, for degree = 2 or larger, nodes <-- group2
0058 % please note: p1*mdeg1 + p2*mdeg2 = mdeg
0059 % so mdeg2 = (mdeg-p1*mdeg1)/p2.
0060 
0061 % calculate the shift factors
0062 p1 = 0.1; mdeg1 = 1.0;
0063 mdeg2 = (mdeg-p1*mdeg1)/(1-p1); % the average degree for group2
0064 
0065 ks = ones(1,N)*mdeg1;
0066 N1 = round(p1*N); N2 = N-N1;  % nodes count
0067 % please note: for matlab Geometric distribution
0068 % prob(x=k) = (1-p)^k*p; k = 0,1,2,....
0069 % E(X) = 1/p -1;
0070 % here, we use this to generate X of {0,1,2,...} then shift it to
0071 % Y = X+2, so that Y of {2,3,4,...}
0072 % E(Y) = 1/p -1 +2 = 1/p+1; on the other side, E(y) = mdeg2;
0073 % therefore p has to be set as: 1/p = mdeg2-1.
0074 ks(N1+1:N) = geornd( 1/(mdeg2 - 1), [1 N2] )+2; % previously set mdeg2-2, it is wrong.
0075 yout = ks(randperm(N));