%computes indices of each adjacency matrix and numbers of neighbors for each block.
%input: msa number and neighbor distance; returns ai, aj vectors of indices of 
%same race neighbors for each adjacency matrix as well as values and sizes

%See Detailed Description (word document) for more details about each of
%the steps labeled below

function [aij1,aij2,aij3,aij4,neigh,sizes,blockid] = neighbors(msa,distance)
global avgdeg avgdegsr

warning off MATLAB:dividebyzero

%Step One: Preliminaries
msafile='msa_number.txt';
msafile=strrep(msafile,'number',msa);
[id,lat,lng,race]=textread(msafile,'%s %f %f %u','delimiter',',');
id=strrep(id,'"','');
id=str2num(cell2mat(id));

k=ones(4,1);
numneigh=zeros(4,length(lat));
blockid=zeros(4,length(id));
aij1=[];
aij2=[];
aij3=[];
aij4=[];
block=zeros(2,4);
%Step Two: Neighbor calculations; iterates through blocks, calculates
%distances and counts both "neighbor" and "same-race neighbor"
%relationships between blocks
for i=1:length(lat)
    block(1,race(i))=block(1,race(i))+1;
    blockid(race(i),block(1,race(i)))=id(i);
    block(2,:)=block(1,:);
    for j=(i+1):length(lat)
        block(2,race(j))=block(2,race(j))+1;
        if abs(lat(i)-lat(j)) <= .03*distance && abs(lng(i)-lng(j)) <= .03*distance
            d=6378.1*acos(sin(2*pi*lat(i)/360)*sin(2*pi*lat(j)/360)+cos(2*pi*(lng(i)-lng(j))/360)*cos(2*pi*lat(i)/360)*cos(2*pi*lat(j)/360));
            if d <= distance
                numneigh(race(i),block(1,race(i)))=numneigh(race(i),block(1,race(i)))+1;
                numneigh(race(j),block(2,race(j)))=numneigh(race(j),block(2,race(j)))+1;
                if race(i)==race(j) 
                    if race(i)==1
                        aij1(:,k(1))=block(:,1);
                        k(1)=k(1)+1;
                    elseif race(i)==2
                        aij2(:,k(2))=block(:,2);
                        k(2)=k(2)+1;
                    elseif race(i)==3
                        aij3(:,k(3))=block(:,3);
                        k(3)=k(3)+1;
                    else
                        aij4(:,k(4))=block(:,4);
                        k(4)=k(4)+1;
                    end
                end
            end
        end
    end
end

sizes=block(1,:);       %number of blocks of each race (1 x 4)

% shrinking sizes to exclude zeroes
zeroneigh(1)=size(find(numneigh(1,1:sizes(1))==0),2);
zeroneigh(2)=size(find(numneigh(2,1:sizes(2))==0),2);
zeroneigh(3)=size(find(numneigh(3,1:sizes(3))==0),2);
zeroneigh(4)=size(find(numneigh(4,1:sizes(4))==0),2);
sizes=sizes-zeroneigh;

% save average degree for all races and for each race
degrees = sum(numneigh,2);
avgdeg(1) = degrees(1) / sizes(1);
avgdeg(2) = degrees(2) / sizes(2);
avgdeg(3) = degrees(3) / sizes(3);
avgdeg(4) = degrees(4) / sizes(4);

numneigh=1./numneigh;   % infinities left in 
% flip the ones
aijtemp=flipud(aij1);
aij1=cat(2,aij1,aijtemp);
aijtemp=flipud(aij2);
aij2=cat(2,aij2,aijtemp);
aijtemp=flipud(aij3);
aij3=cat(2,aij3,aijtemp);
aijtemp=flipud(aij4);
aij4=cat(2,aij4,aijtemp);
%match neighbor numbers with ai&aj's
for j=1:size(aij1,2)
    neigh(1,j)=numneigh(1,aij1(1,j));
end
for j=1:size(aij2,2)
    neigh(2,j)=numneigh(2,aij2(1,j));
end
for j=1:size(aij3,2)
    neigh(3,j)=numneigh(3,aij3(1,j));
end
for j=1:size(aij4,2)
    neigh(4,j)=numneigh(4,aij4(1,j));
end

avgdegsr(1) = length(aij1) / sizes(1);
avgdegsr(2) = length(aij2) / sizes(2);
avgdegsr(3) = length(aij3) / sizes(3);
avgdegsr(4) = length(aij4) / sizes(4);