function [tie] = establish_ties_gss(h, rie)

% function [tie] = establish_ties_gss(h, rie)
%
% Simulates a network of relations with homophily. Inputs are the homophily 
% coefficent (h), a matrix of race, income, education and no of ties(rie). 
% Output is a tie matrix. Target number of ties for a person is given by 
% 'numprobs' from GSS, stored in the 8th column of rie. Homophily is based 
% on a social distance metric (i.e., when homophily occurs, a random tie is 
% established with an individual j that is among the n people closest to the 
% individual i in terms of the social distance metric.) n is arbitrarily
% chosen to be 20*target # of relations for an individual (so that in the case
% of perfect homophily, the individual has sufficient number of ingroup 
% members to choose from. 
% The weigths of race, income and education in the social distance equation
% are estimated using the homophily estimates in McPherson, Smith-Lovin and
% Brashears (2006).
%
% INPUT:    h    - Homophily coefficient
%           rie  - Node #, Race, Income, Education, Income Group ID,
%                  Education Group ID Matrix, Number of Contacts, Number of 
%                  friends to discuss personal or other problems
% OUTPUT:   tie  - The tie matrix (0/1)
%           
% Paul DiMaggio and Filiz Garip (c)
% Department of Sociology
% Princeton University, Harvard University
% Created: May 27, 2008
% File Name: establish_ties_gss.m


ps    = size(rie,1);
tgt   = rie(:,8);


% CREATE A SOCIAL DISTANCE MATRIX AND DETERMINE IN-GROUP MEMBERS

socd         = zeros(ps,ps);     % Social Distance Matrix
ingr         = zeros(ps,ps);     % Ingroup Matrix                            

% Create a new matrix where race, income and educ are standardized into [0,1] range
ries        = rie(:,:);          
ries(:,3)   = (ries(:,3) - min(ries(:,3))) / range(ries(:,3));
ries(:,4)   = (ries(:,4) - min(ries(:,4))) / range(ries(:,4));

% Social distance weights based on McPherson et al. (2006), Table 3, p.360
% They use GSS 2004 Networks Data (NUMGIVEN variable) to measure the
% homophily in social networks. We compute the weights as follows:
% w = 1-(Heterogeneity in networks/Heterogeneity in population)
% As there is no information on income homophily, we assume that the weight
% of income is the same as that of education in social distance
% computation.

rw = 1 - (0.09/0.53);
ew = 1 - (1.48/3.17);
iw = ew;

% Compute social distance among each pair (i,j)
for i = 1:ps
    for j = 1:ps
        socd(i,j)  = sqrt((rw*(ries(i,2)-ries(j,2)))^2 + (iw*(ries(i,3)-ries(j,3)))^2 + (ew*(ries(i,4)-ries(j,4)))^2);
    end
end

% Assign (20*target #) of in-group members to each individual i based on 
% social distance. Note that the number of in-group members is selected 
% arbitrarily to be twenty times the target number of ties for that person. 
% This is done to account for the possibility that some ties will be infeasible 
% given a maximum number of ties allowed for each individual. Also note that 
% because social distance to self is zero, self is initially selected as 
% an ingroup member. The code then corrects for this, and sets ingr(i,i)=0 for all i.

for i = 1:ps
    [Y, I] = sort(socd(:,i));           % Sort social distance matrix (socd) by column i
    for j = 1:(20*tgt(i,1))
        ingr(I(j,1),i) = 1;             % Select tgt(i,1) individuals from the sorted socd, and flag them as in-group members 
    end
end
ingr = ingr - eye(ps);                  % Exclude self from ingroup members


% ESTABLISH NETWORK TIES

tie         = eye(ps);              % Matrix of ties b/w nodes (0/1) - tie to oneself is set to 1.
c_h         = 0;                    % Count homophily occurence
ct          = zeros(ps,1);          % Vector of current ties per person

avn         = zeros(ps,2);          % Matrix of available nodes: Col (1)0/1 for availability, 
                                    % (2) Node number of available nodes
avn(:,1)    = (ct<tgt);             % Mark the available nodes - initially all available
avn(:,2)    = [1:ps]';              % Give them numbers

noingr      = zeros(ps,1);          % A vector that flags nodes with no more potential in-group ties. 
avnh        = avn;                  % Available nodes under homophily (i.e., nodes with potential in-group ties available)

while (sum(avn(:,1))>0)             % While there are available nodes
   
    % If homophily bias occurs, select a random node from in-group.
    
    if (rand<=h)             
        if (sum(avnh(:,1))== 0)
            break;
        end
        u       = unidrnd(sum(avnh(:,1)));                  % Pick a random node from available nodes under homophily
        Iu      = find(avnh(:,2) == u);                     % Find that node
        Ig      = find(ingr(:,Iu) == 1 & avn(:,1) == 1);    % Index all the group members of the selected node that are available 
                                                            % (i.e., nodes that have not yet reached their target number of relations) 
                                                          
        % Establish a tie with an ingroup member
        
        target     = tie(Iu,Ig);                    % Submatrix of the selected node's ties to in-group members
        no_tie_ind = find(target == 0);             % Index the potential ties (i.e., not yet established)
        no_tie     = size(Ig,1) - sum(sum(target)); % Count the number of potential ties
        if (no_tie ~= 0)                                        % If number of potential ties is not zero
            t                               = unidrnd(no_tie);  % Pick a random node number
            tie(Iu,Ig(no_tie_ind(t)))       = 1;                % Establish a tie
            tie(Ig(no_tie_ind(t)),Iu)       = 1;
            ct(Iu)                          = ct(Iu) + 1;       % Increment number of current ties by one for both nodes
            if (ct(Iu) > tgt(Iu))
                display('Exceeded target tie number');
            end
            ct(Ig(no_tie_ind(t)))           = ct(Ig(no_tie_ind(t))) + 1;
            if (ct(Ig(no_tie_ind(t))) > tgt(Ig(no_tie_ind(t))))
                display('Exceeded target tie number');
            end;
            c_h                             = c_h + 1;          % Increment count of ties with homophily
        else 
            % disp('No more in-group ties available');               
            noingr(Iu,1) = 1;       % Flag this node to be dropped from the avaliable node list
        end
     
    % If homophily does not occur
    
    else
        u          = unidrnd(sum(avn(:,1)));        % Pick a random node from all available nodes 
        Iu         = find(avn(:,2) == u);           % Find that node
        It         = find(avn(:,1)==1);             % Index all available nodes
        target     = tie(Iu,It);                    % Submatrix of the selected node's ties to all available members
        no_tie_ind = find(target == 0);             % Index the potential ties 
        no_tie     = size(It,1) - sum(sum(target)); % Count the number of potential ties
        
        if (no_tie ~= 0)                                        % If number of potential ties is not zero
            t                           = unidrnd(no_tie);      % Pick a random node number
            tie(Iu,It(no_tie_ind(t)))       = 1;                    % Establish a tie
            tie(It(no_tie_ind(t)),Iu)       = 1;
            ct(Iu)                      = ct(Iu) + 1;           % Increment number of current ties by one for both nodes
            ct(It(no_tie_ind(t)))           = ct(It(no_tie_ind(t))) + 1;
            if (ct(Iu) > tgt(Iu))
                display('Exceeded target tie number');
            end
            if (ct(It(no_tie_ind(t))) > tgt(It(no_tie_ind(t))))
                display('Exceeded target tie number');
            end;
            
        else
            % disp('No more nodes available');
            break;                     % Break if there are no nodes available
        end
    end
    avn         = zeros(ps,2);          % Clear prior values for available nodes
    avn(:,1)    = (ct<tgt);             % Update the available nodes
    [I J]       = find(avn(:,1) == 1);  % Index the available nodes
    avn(I,2)    = [1:size(I,1)]';       % Give node numbers
    
    avnh        = zeros(ps,2);          
    
    avnh(:,1)   = avn(:,1)- noingr + (avn(:,1) < noingr);     % Update the available nodes under homophily 
                                                              % (i.e. discard nodes with no in-group ties available)
                                                              % Avoid getting negative values once avn(:,1)< noingr
    
    [I J]       = find(avnh(:,1) == 1); % Index the available nodes
    avnh(I,2)   = [1:size(I,1)]';       % Give node numbers
    
end

disp('Ties that could not be established');     % Typically due to homophily constraints
disp(sum(tgt-ct));
 
tie = tie - eye(ps);                            % Set tie to oneself to zero.


return;