program centpow
version 11
//Revised 12jan2011
syntax anything [, BONacich NORMalize beta(real 0) SYMmetrize(str) saveas(str)]
quietly {
preserve

//Check that beta is positive
if `beta' < 0 {
	noisily: display "{bf:ERROR -} The {it:beta} parameter must be positive."
	noisily: display " "
	}
error (`beta' < 0)

//Import data
clear
insheet using `anything'
mata: network = st_data(.,.)
mata: _diag(network, 0)

//Check that matrix is square
mata: row = rows(network)
mata: col = cols(network)
mata: st_numscalar("r(row)", row)
mata: st_numscalar("r(col)", col)
if r(row) ~= r(col) {
	noisily: display "{bf:NOTICE -} The matrix is not square."
	error (r(row) ~= r(col))
	}

//Symmetrize
mata: symmetric = issymmetric(network)
mata: st_numscalar("r(symmetric)", symmetric)
if r(symmetric) == 1 & "`symmetrize'" ~= "" {
	noisily: display "{bf:NOTICE -} The matrix is already symmetric."
	noisily: display "         It will be analyzed in its present form."
	noisily: display " "
	}
if r(symmetric) == 0 {
	if "`symmetrize'" == "upper" {
		mata: _uppertriangle(network)
		mata: network = network'
		mata: _makesymmetric(network)
		}
	else if "`symmetrize'" == "lower" {
		mata: _lowertriangle(network)
		mata: _makesymmetric(network)
		}
	else if "`symmetrize'" == "sum" | "`symmetrize'" == "" {
		mata: upper = uppertriangle(network)
		mata: upper = upper'
		mata: _makesymmetric(upper)
		mata: lower = lowertriangle(network)
		mata: _makesymmetric(lower)
		mata: network = upper + lower
		}
	else {
		noisily: display "{bf:NOTICE -} The {it:symmetrize} option is incorrectly specified."
		error (r(symmetric) == 0)
		}
	}

//If normalization is requested, check that matrix is binary
if "`normalize'" ~= "" {
	mata: max = max(network)
	mata: st_numscalar("r(max)", max)
	if r(max) > 1 {
		noisily: display "{bf:ERROR -} Normalized values can only be computed for binary matrices."
		noisily: display " "
		}
	error (r(max) > 1)
	}

//Compute indices
mata: degree = rowsum(network)
mata: altercent = network[,] * degree
mata: invdeg = 1 :/ degree
mata: _editmissing(invdeg,0)
mata: alterpow = network[,] * invdeg
if "`bonacich'" ~= "" {
	mata: y = I(rows(network))
	mata: z = J(rows(network),1,1)
	mata: symeigensystem(network,eigenvector=.,eigenvalue=.)
	mata: dimension = eigenvalue[1,1]/eigenvalue[1,2]
	mata: maxbeta = 1 / eigenvalue[1,1]
	if `beta' == 0 mata: beta = Re(maxbeta * .99)
		else mata: beta = `beta'
	mata: betacent = (luinv(y - (beta * network)))*network*z
	mata: betapow = (luinv(y - (-beta * network)))*network*z
	mata: eigencent = abs(eigenvector[.,1])
	
	//Flag violations of beta centrality/power assumptions
	mata: st_numscalar("r(beta)", beta)
	mata: st_numscalar("r(maxbeta)", Re(maxbeta))
	mata: st_numscalar("r(dimension)", Re(dimension))

	mata: x = network			//Flag networks with multiple components, r(minpath) = 0
	mata: y = network
	mata: z = network
	local path = 1
	while `path' < (5 - 1) {
		mata: y = y * x
		mata: z = z + y
		local path = `path' + 1
		}
	mata: st_numscalar("r(minpath)", min(z))

	if "`bonacich'" ~= "" noisily: display "Beta parameter: " r(beta)

	if r(beta) > r(maxbeta) {
		noisily: display"{bf:WARNING -} Beta centrality and power should be interpreted with"
		noisily: display "caution because the selected value of {it:beta} (" r(beta) ") exceeds"
		noisily: display "the maximum allowable value of " round(r(maxbeta), .001) "."
		noisily: display " "
		}
	if r(beta) < .05 {
		noisily: display "{bf:WARNING -} Beta centrality and power results are nearly the same as ordinary"
		noisily: display "degree centrality because the value of {it:beta} (" r(beta) ") is nearly zero."
		noisily: display " "
		}
	if r(dimension) < 2 {
		noisily: display "{bf:WARNING -} Eigenvector centrality, beta centrality and beta power should be"
		noisily: display "interpreted with caution because the network's largest eigenvalue is only"
		noisily: display round(r(dimension), .01) " times the size of the 2nd largest eigenvalue.""
		noisily: display " "
		}
	if r(minpath) == 0 {
		noisily: display "{bf:WARNING -} Eigenvector centrality, beta centrality and beta power should be"
		noisily: display "interpreted with caution because the network contains multiple components."
		noisily: display " "
		}
	}
if "`normalize'" ~= "" {
	mata: size = rows(network)
	mata: cent_max = (size - 1) * (size - 1)
	mata: cent_min = 2
	mata: n_altercent = (altercent :- cent_min) :/ (cent_max - cent_min)

	//Correct normalization for isolated dyads
	mata: neg = min(n_altercent)
	mata: st_numscalar("r(neg)", neg)
	if r(neg) < 0 {
		mata: n_altercent = (altercent :- 1) :/ (cent_max - 1)
		}

	mata: pow_max = size - 1
	mata: pow_min = 1 / (size - 1)
	mata: n_alterpow = (alterpow :- pow_min) :/ (pow_max - pow_min)
	}

//Post results to Stata
clear
if "`bonacich'" ~= "" & "`normalize'" ~= "" mata: results = (degree, altercent, alterpow, n_altercent, n_alterpow, betacent, betapow, eigencent)
else if "`bonacich'" ~= "" & "`normalize'" == "" mata: results = (degree, altercent, alterpow, betacent, betapow, eigencent)
else if "`bonacich'" == "" & "`normalize'" ~= "" mata: results = (degree, altercent, alterpow, n_altercent, n_alterpow)
else if "`bonacich'" == "" & "`normalize'" == "" mata: results = (degree, altercent, alterpow)
mata: st_matrix("results", results)
svmat results
rename results1 degree
rename results2 altercent
rename results3 alterpow
if "`bonacich'" ~= "" & "`normalize'" ~= "" rename results4 n_altercent
if "`bonacich'" ~= "" & "`normalize'" ~= "" rename results5 n_alterpow
if "`bonacich'" ~= "" & "`normalize'" ~= "" rename results6 betacent
if "`bonacich'" ~= "" & "`normalize'" ~= "" rename results7 betapow
if "`bonacich'" ~= "" & "`normalize'" ~= "" rename results8 eigencent
if "`bonacich'" ~= "" & "`normalize'" == "" rename results4 betacent
if "`bonacich'" ~= "" & "`normalize'" == "" rename results5 betapow
if "`bonacich'" ~= "" & "`normalize'" == "" rename results6 eigencent
if "`bonacich'" == "" & "`normalize'" ~= "" rename results4 n_altercent
if "`bonacich'" == "" & "`normalize'" ~= "" rename results5 n_alterpow

//Save results and end
if "`saveas'" ~= "" save "`saveas'"
else save centpow
restore
}
end