igraphHomophily.R

#####
### Except where otherwise indicated, the code below is written by Benjamin Lind under CC: BY-NC-SA 3.0
### If you use or develop this code, please provide proper attribution.
### You may contact me at lind.benjamin//at//gmail//d0t//com
### You may also find me @benjaminlind on Twitter
#####

### To load these functions into R, enter:
###     source("http://pastebin.com/raw.php?i=s2ZcqxSm")

### Functions included in this file:
###     Krackhardt's E-I index
###     Graph generation based on homophily
###         Homophily graphs can be created according to general tie probabilities within and between groups.
###         Homophily graphs can also created according to tie probabilities within and between specific groups.
###         These functions can be used for any categorical attribute, including blocks or communities.
###     Tie probability calculators within and between groups.
###         The values produced from these functions can be used in the homophily graph generators.
###     Gould and Fernandez's (1989) brokerage measurements
###         Offers corrections to brokerage() on undirected graphs in library(sna)
###         Includes an attribute permutation test
###     Additional functions not related to homophily, but otherwise useful and not in igraph
###         Weak transitivity as measured for a directed network
###         UMAN random graph generators. They can make random graphs conditioned on the dyad census.
###         A dyad rewired graph. Like igraph::rewire(), but controls for reciprocity.
###         Vazquez random graph generator.
###         Krackhardt's (1994) measurements of hierarchical organizations. These include:
###             Hierarchy
###             Efficiency
###             Connectedness
###             LUBness.
###                 The code for LUBness below was initially written by Butts (2014) for library(sna).
###                 I have adapted the code here for igraph and Rcpp.

library(Rcpp)
library(igraph)
library(compiler)

cppFunction('
    double ei_index_cpp(IntegerMatrix el, IntegerVector attrib){
        double incount = 0;
        int n = el.nrow();
        for (int i = 0; i < n; ++i){
            int attrib1 = attrib[el(i, 0) - 1];
            int attrib2 = attrib[el(i, 1) - 1];
            if (attrib1 == attrib2){
                incount += 1;
                }
            }
        double excount = n - incount;
        double out = (excount - incount) / n;
        return out;
        }
')

cppFunction('
    LogicalMatrix homophily_game_cpp_general(IntegerVector m, NumericVector p, bool dir, bool loop){
        int n = m.size();
        LogicalMatrix out(n, n);
        bool samegroup = false;
        double tp = 0.0;
        int addedgeint = 0;
        double addedgedoub = 0.0;
        bool addedge = false;

        if (dir & (loop == false)){
            for (int i = 0; i < n; ++i){
                for (int j = 0; j < n; ++j){
                    if (i != j){
                        samegroup = (m[i] == m[j]);
                        if (samegroup){
                            tp = p[0];
                            }
                        else{
                            tp = p[1];
                            }
                        addedgeint = rand() % 10000000 + 1;
                        addedgedoub = addedgeint / 10000000.0;
                        addedge = addedgedoub < tp;
                        if (addedge){
                            out(i, j) = true;
                            }
                        else{
                            out(i, j) = false;
                            }
                        }
                    else{
                        out(i, j) = false;
                        }
                    }
                }
            }
        else if (dir & loop){
            for (int i = 0; i < n; ++i){
                for (int j = 0; j < n; ++j){
                    samegroup = (m[i] == m[j]);
                    if (samegroup){
                        tp = p[0];
                        }
                    else{
                        tp = p[1];
                        }
                    addedgeint = rand() % 10000000 + 1;
                    addedgedoub = addedgeint / 10000000.0;
                    addedge = addedgedoub < tp;
                    if (addedge){
                        out(i, j) = true;
                        }
                    else{
                        out(i, j) = false;
                        }
                    }
                }
            }
        else if ((dir == false) & (loop == false)){
            for (int i = 0; i < n; ++i){
                for (int j = i; j < n; ++j){
                    if (i != j){
                        samegroup = (m[i] == m[j]);
                        if (samegroup){
                            tp = p[0];
                            }
                        else{
                            tp = p[1];
                            }
                        addedgeint = rand() % 10000000 + 1;
                        addedgedoub = addedgeint / 10000000.0;
                        addedge = addedgedoub < tp;
                        if (addedge){
                            out(i, j) = true;
                            out(j, i) = true;
                            }
                        else{
                            out(i, j) = false;
                            out(j, i) = false;
                            }
                        }
                    else{
                        out(i, j) = false;
                        }
                    }
                }
            }
        else{
            for (int i = 0; i < n; ++i){
                for (int j = i; j < n; ++j){
                    samegroup = (m[i] == m[j]);
                    if (samegroup){
                        tp = p[0];
                        }
                    else{
                        tp = p[1];
                        }
                    addedgeint = rand() % 10000000 + 1;
                    addedgedoub = addedgeint / 10000000.0;
                    addedge = addedgedoub < tp;
                    if (addedge){
                        out(i, j) = true;
                        out(j, i) = true;
                        }
                    else{
                        out(i, j) = false;
                        out(j, i) = false;
                        }
                    }
                }
            }
        return out;
        }
')

# To make a faster version of the next function, use the probabilities to estimate the number of edges in each pairing.
# Then, for each pairing, randomly select a pair, verify don't already exist in the edge list, then add to the edge list until the quota is met.
# (A similar method could be adopted that integrates dyads into each group pairing.)

cppFunction('
    LogicalMatrix homophily_game_cpp_mixed(IntegerVector m, NumericMatrix p, bool dir, bool loop){
        int n = m.size();
        LogicalMatrix out(n, n);
        double tp = 0.0;
        int addedgeint = 0;
        double addedgedoub = 0.0;
        bool addedge = false;

        if (dir & (loop == false)){
            for (int i = 0; i < n; ++i){
                for (int j = 0; j < n; ++j){
                    if (i != j){
                        tp = p(m[i] - 1, m[j] - 1);
                        addedgeint = rand() % 10000000 + 1;
                        addedgedoub = addedgeint / 10000000.0;
                        addedge = addedgedoub < tp;
                        if (addedge){
                            out(i, j) = true;
                            }
                        else{
                            out(i, j) = false;
                            }
                        }
                    else{
                        out(i, j) = false;
                        }
                    }
                }
            }
        else if (dir & loop){
            for (int i = 0; i < n; ++i){
                for (int j = 0; j < n; ++j){
                    tp = p(m[i] - 1, m[j] - 1);
                    addedgeint = rand() % 10000000 + 1;
                    addedgedoub = addedgeint / 10000000.0;
                    addedge = addedgedoub < tp;
                    if (addedge){
                        out(i, j) = true;
                        }
                    else{
                        out(i, j) = false;
                        }
                    }
                }
            }
        else if ((dir == false) & (loop == false)){
            for (int i = 0; i < n; ++i){
                for (int j = i; j < n; ++j){
                    if (i != j){
                        tp = p(m[i] - 1, m[j] - 1);
                        addedgeint = rand() % 10000000 + 1;
                        addedgedoub = addedgeint / 10000000.0;
                        addedge = addedgedoub < tp;
                        if (addedge){
                            out(i, j) = true;
                            out(j, i) = true;
                            }
                        else{
                            out(i, j) = false;
                            out(j, i) = false;
                            }
                        }
                    else{
                        out(i, j) = false;
                        }
                    }
                }
            }
        else{
            for (int i = 0; i < n; ++i){
                for (int j = i; j < n; ++j){
                    tp = p(m[i] - 1, m[j] - 1);
                    addedgeint = rand() % 10000000 + 1;
                    addedgedoub = addedgeint / 10000000.0;
                    addedge = addedgedoub < tp;
                    if (addedge){
                        out(i, j) = true;
                        out(j, i) = true;
                        }
                    else{
                        out(i, j) = false;
                        out(j, i) = false;
                        }
                    }
                }
            }
        return out;
        }
')

cppFunction('
    NumericVector group_density_match_cpp(IntegerMatrix el, IntegerVector attrib, bool dir, bool loop){
        double incount = 0;
        int n = el.nrow();
        int attrib1 = 0;
        int attrib2 = 0;

        for (int i = 0; i < n; ++i){
            attrib1 = attrib[el(i, 0) - 1];
            attrib2 = attrib[el(i, 1) - 1];
            if (attrib1 == attrib2){
                incount += 1;
                }
            }
        double excount = n - incount;
        n = attrib.size();
        double nposse = 0;
        if (dir & loop){
            nposse = pow(n, 2);
            }
        else if (dir & (loop == false)){
            nposse = n * (n - 1);
            }
        else if ((dir == false) & loop){
            nposse = ((n * (n - 1)) / 2) + n;
            }
        else{
            nposse = (n * (n - 1)) / 2;
            }

        int mmax = 0;
        for (int i = 0; i < n; ++i){
            if (attrib[i] > mmax){
                mmax = attrib[i];
                }
            }

        IntegerVector gsize(mmax);
        for (int i = 0; i < mmax; ++i){
            gsize[i] = 0;
            }
        for (int i = 0; i < n; ++i){
            gsize[attrib[i] - 1] += 1;
            }

        double npossein = 0.0;
        if (dir & loop){
            for (int i = 0; i < mmax; ++i){
                npossein += pow(gsize[i], 2);
                }
            }
        else if (dir & (loop == false)){
            for (int i = 0; i < mmax; ++i){
                npossein += gsize[i] * (gsize[i] - 1);
                }
            }
        else if ((dir == false) & loop){
            for (int i = 0; i < mmax; ++i){
                npossein += ((gsize[i] * (gsize[i] - 1)) / 2) + gsize[i];
                }
            }
        else{
            for (int i = 0; i < mmax; ++i){
                npossein += (gsize[i] * (gsize[i] - 1)) / 2;
                }
            }

        NumericVector out(2);
        out[0] = incount / npossein;
        out[1] = excount / (nposse - npossein);
        out.names() = CharacterVector::create("internal", "external");
        return out;
        }
')

cppFunction('
    NumericMatrix group_density_mixture_cpp(IntegerMatrix el, IntegerVector attrib, bool dir, bool loop){
        int mmax = 0;
        int n = attrib.size();

        for (int i = 0; i < n; ++i){
            if (attrib[i] > mmax){
                mmax = attrib[i];
                }
            }

        NumericMatrix ecounts(mmax, mmax);
        int attrib1 = 0;
        int attrib2 = 0;
        int ne = el.nrow();

        for (int i = 0; i < mmax; ++i){
            for (int j = i; j < mmax; ++j){
                ecounts(i, j) = 0.0;
                ecounts(j, i) = 0.0;
                }
            }

        if (dir){
            for (int i = 0; i < ne; ++i){
                attrib1 = attrib[el(i, 0) - 1] - 1;
                attrib2 = attrib[el(i, 1) - 1] - 1;
                ecounts(attrib1, attrib2) += 1.0;
                }
            }
        else{
            for (int i = 0; i < ne; ++i){
                attrib1 = attrib[el(i, 0) - 1] - 1;
                attrib2 = attrib[el(i, 1) - 1] - 1;
                ecounts(attrib1, attrib2) += 1;
                if (attrib1 != attrib2){
                    ecounts(attrib2, attrib1) += 1.0;
                    }
                }
            }

        NumericVector gsize(mmax);
        for (int i = 0; i < mmax; ++i){
            gsize[i] = 0.0;
            }
        for (int i = 0; i < n; ++i){
            gsize[attrib[i] - 1] += 1;
            }
        if (dir & loop){
            for (int i = 0; i < mmax; ++i){
                for (int j = i; j < mmax; ++j){
                    if (i != j){
                        ecounts(i, j) = ecounts(i, j) / (gsize[i] * gsize[j]);
                        ecounts(j, i) = ecounts(j, i) / (gsize[j] * gsize[i]);
                        }
                    else{
                        ecounts(i, j) = ecounts(i, j) / pow(gsize[i], 2);
                        }
                    }
                }
            }
        else if (dir & (loop == false)){
            for (int i = 0; i < mmax; ++i){
                for (int j = i; j < mmax; ++j){
                    if (i != j){
                        ecounts(i, j) = ecounts(i, j) / (gsize[i] * gsize[j]);
                        ecounts(j, i) = ecounts(j, i) / (gsize[j] * gsize[i]);
                        }
                    else{
                        ecounts(i, j) = ecounts(i, j) / (gsize[i] * (gsize[i] - 1));
                        }
                    }
                }
            }
        else if ((dir == false) & loop){
            for (int i = 0; i < mmax; ++i){
                for (int j = i; j < mmax; ++j){
                    if (i != j){
                        ecounts(i, j) = ecounts(i, j) / (gsize[i] * gsize[j]);
                        ecounts(j, i) = ecounts(j, i) / (gsize[j] * gsize[i]);
                        }
                    else{
                        ecounts(i, j) = ecounts(i, j) / (((gsize[i] * (gsize[i] - 1)) / 2) + gsize[i]);
                        }
                    }
                }
            }
        else{
            for (int i = 0; i < mmax; ++i){
                for (int j = i; j < mmax; ++j){
                    if (i != j){
                        ecounts(i, j) = ecounts(i, j) / (gsize[i] * gsize[j]);
                        ecounts(j, i) = ecounts(j, i) / (gsize[j] * gsize[i]);
                        }
                    else{
                        ecounts(i, j) = ecounts(i, j) / ((gsize[i] * (gsize[i] - 1)) / 2);
                        }
                    }
                }
            }

        return ecounts;
        }
')

cppFunction('
    IntegerMatrix replace_multiple_ties_cpp(IntegerMatrix el, IntegerVector sinds, IntegerVector rinds, bool dir, bool loop){
        int nrow = el.nrow();
        int ncol = el.ncol();
        IntegerMatrix out(nrow, ncol);
        out(0, 0) = el(0, 0);
        out(0, 1) = el(0, 1);
        for (int i = 1; i < nrow; ++i){
            bool reprow = false;
            if ((loop == false) & (el(i, 0) == el(i, 1))){
                reprow = true;
                }
            int j = 0;
            while ((j < i) & (reprow == false)){
                if ((el(i, 0) == out(j, 0)) & (el(i, 1) == out(j, 1))){
                    reprow = true;
                    }
                if ((dir == false) & ((el(i, 0) == out(j, 1)) & (el(i, 1) == out(j, 0)))){
                    reprow = true;
                    }
                j += 1;
                }
            out(i, 0) = el(i, 0);
            out(i, 1) = el(i, 1);
            while (reprow == true){
                out(i, 0) = sinds[rand() % sinds.size()];
                out(i, 1) = rinds[rand() % rinds.size()];
                bool newreprow = false;
                if ((loop == false) & (out(i, 0) == out(i, 1))){
                    newreprow = true;
                    }
                int j = 0;
                while ((j < i) & (newreprow == false)){
                    if ((out(i, 0) == out(j, 0)) & (out(i, 1) == out(j, 1))){
                        newreprow = true;
                        }
                    if ((dir == false) & ((out(i, 0) == out(j, 0)) & (out(i, 1) == out(j, 1)))){
                        newreprow = true;
                        }
                    j += 1;
                    }
                reprow = newreprow;
                }
            }
        return(out);
        }
')

homophily_cpp <- list(ei = ei_index_cpp, group_density = list(general = group_density_match_cpp, mixed = group_density_mixture_cpp), games = list(general = homophily_game_cpp_general, mixed = homophily_game_cpp_mixed), assist = list(replace_multiple_ties = replace_multiple_ties_cpp))

rm(ei_index_cpp, group_density_match_cpp, group_density_mixture_cpp, homophily_game_cpp_general, homophily_game_cpp_mixed, replace_multiple_ties_cpp)

ei_index <- function(g, value){
    # _g_ is an igraph 'graph' object
    # _value_ is a categorical vertex attribute
    if (any(is.na(value)))
        stop("Missing values are not allowed for _value_.")
    if (vcount(g) != length(value))
        stop("The length of _value_ does not equal the number of vertices in _g_.")
    g.el <- get.edgelist(g, names = FALSE)
    allvalues <- unique(value)
    value.int <- match(value, allvalues)
    return(homophily_cpp$ei(g.el, value.int))
    }
ei_index <- cmpfun(ei_index)

homophily_game <- function(groups.n, tie.prob, directed = TRUE, loops = FALSE, sparse = FALSE){
    # _groups.n_ is a vector of integers indicating the size of each group.
    # _tie.prob_ needs to be either a matrix or a vector of length two.
        # If a vector
            # The first value needs to refer to the tie probability within groups
            # The second value needs to refer to the tie probability between groups
        # If a matrix
            # Values are proportional to the tie probability within and between groups.
            # It must be square
        # Note that the tie probability is for the pairings. The denominator varies depending on group size.
    # _directed_ indicates whether or not the simulated network should be directed.
    # _loops_ indicates whether or not loops should be allowed.
    # Set _sparse_ equal to TRUE for large (> ~1000 vertices) and sparse networks ( < ~0.01 density).
        # The default argument (FALSE) create matrices that occupy considerable memory for large networks.
            # It works great for small networks, as indexing is very fast.
        # The sparse algorithm involves creating an edge list with randomly selected nodes.
            # This technique introduces the possibility for multiple edges to occur.
            # To correct the possibility of multiple edges, the algorithm loops through the edge list to replace them.

    groups.n <- round(groups.n)
    s <- length(groups.n)
    if (s < 1)
        stop("There must be at least two groups indicated by _groups.n_.")
    if (any(tie.prob < 0) | any(tie.prob > 1))
        stop("_tie.prob_ must have all values be greater than zero and less than one to indicate the desnity within and between groups.")
    if (any(is.na(tie.prob)))
        stop("Missing data in _tie.prob_ is not supported.")
    if ((length(tie.prob) != 2) & (length(tie.prob) != s^2))
        stop("_tie.prob_ needs either two values or a square density matrix.")
    if (is.logical(directed) == FALSE)
        stop("_directed_ must be a logical argument to indicate directionality of the desired graph.")

    n <- sum(groups.n)
    group.assignments <- unlist(sapply(1:s, function(x) rep(x, times = groups.n[x])))
    if (length(tie.prob) > 2){
        if (nrow(tie.prob) != ncol(tie.prob))
            stop("Number of columns does not equal the number of rows provided in _tie.prob_")
        if (nrow(tie.prob) != s)
            stop("Number of rows and columns does not equal the number of groups provided in _tie.prob_.")
        if (sparse){
            add.inds <- c(0, cumsum(groups.n[1:(s-1)]))
            intergroup.graphs <- lapply(1:s, function(x){
                newg <- erdos.renyi.game(groups.n[x], tie.prob[x, x], type = "gnp", directed = directed, loops = loops)
                newg <- get.edgelist(newg, names = FALSE)
                newg <- as.matrix(newg) + add.inds[x]
                return(newg)
                }
                )
            intergroup.graphs <- do.call(rbind, intergroup.graphs)
            if (directed){
                for (i in 1:s){
                    for (j in 1:s){
                        if (i != j){
                            expties <- round(groups.n[i] * groups.n[j] * tie.prob[i, j])
                            if (expties > 0){
                                send <- add.inds[i] + sample.int(n = groups.n[i], size = expties, replace = TRUE)
                                rec <- add.inds[j] + sample.int(n = groups.n[j], size = expties, replace = TRUE)
                                if (expties == 1)
                                    intergroup.graphs <- rbind(intergroup.graphs, c(send, rec))
                                else{
                                    newg <- matrix(NA, nrow = expties, ncol = 2)
                                    newg[, 1] <- send
                                    newg[, 2] <- rec
                                    possible.mult.ties <- duplicated(send) * duplicated(rec)
                                    if (any(possible.mult.ties)){
                                        iinds <- add.inds[i] + 1:groups.n[i]
                                        jinds <- add.inds[j] + 1:groups.n[j]
                                        newg <- homophily_cpp$assist$replace_multiple_ties(newg, iinds, jinds, directed, loops)
                                        intergroup.graphs <- rbind(intergroup.graphs, newg)
                                        }
                                    else
                                       intergroup.graphs <- rbind(intergroup.graphs, newg)
                                    }
                                }
                            }
                        }
                    }
                }
            else{
                for (i in 1:(s - 1)){
                    for (j in (i + 1):s){
                        expties <- round(groups.n[i] * groups.n[j] * tie.prob[i, j])
                        if (expties > 0){
                            send <- add.inds[i] + sample.int(n = groups.n[i], size = expties, replace = TRUE)
                            rec <- add.inds[j] + sample.int(n = groups.n[j], size = expties, replace = TRUE)
                            if (expties == 1)
                                intergroup.graphs <- rbind(intergroup.graphs, c(send, rec))
                            else{
                                newg <- matrix(nrow = expties, ncol = 2)
                                newg[, 1] <- send
                                newg[, 2] <- rec
                                possible.mult.ties <- duplicated(send) * duplicated(rec)
                                if (any(possible.mult.ties)){
                                    iinds <- add.inds[i] + 1:groups.n[i]
                                    jinds <- add.inds[j] + 1:groups.n[j]
                                    newg <- homophily_cpp$assist$replace_multiple_ties(newg, iinds, jinds, directed, loops)
                                    intergroup.graphs <- rbind(intergroup.graphs, newg)
                                    }
                                else
                                    intergroup.graphs <- rbind(intergroup.graphs, newg)
                                }
                            }
                        }
                    }
                }
            g <- as.data.frame(intergroup.graphs)
            g <- graph.data.frame(g, directed = directed, vertices = data.frame(name = 1:n, group = group.assignments))
            }
        else{
            g <- homophily_cpp$games$mixed(group.assignments, tie.prob, directed, loops)
            g <- graph.adjacency(g, mode = ifelse(directed, "directed", "undirected"), diag = loops)
            }
        }
    else{
        if (sparse){
            add.inds <- c(0, cumsum(groups.n[1:(s-1)]))
            intergroup.graphs <- lapply(1:s, function(x){
                newg <- erdos.renyi.game(groups.n[x], tie.prob[1], type = "gnp", directed = directed, loops = loops)
                newg <- get.edgelist(newg, names = FALSE)
                newg <- as.matrix(newg) + add.inds[x]
                return(newg)
                }
                )
            intergroup.graphs <- do.call(rbind, intergroup.graphs)
            if (directed){
                for (i in 1:s){
                    for (j in 1:s){
                        if (i != j){
                            expties <- round(groups.n[i] * groups.n[j] * tie.prob[2])
                            if (expties > 0){
                                send <- add.inds[i] + sample.int(n = groups.n[i], size = expties, replace = TRUE)
                                rec <- add.inds[j] + sample.int(n = groups.n[j], size = expties, replace = TRUE)
                                if (expties == 1)
                                    intergroup.graphs <- rbind(intergroup.graphs, c(send, rec))
                                else{
                                    newg <- matrix(NA, nrow = expties, ncol = 2)
                                    newg[, 1] <- send
                                    newg[, 2] <- rec
                                    possible.mult.ties <- duplicated(send) * duplicated(rec)
                                    if (any(possible.mult.ties)){
                                        iinds <- add.inds[i] + 1:groups.n[i]
                                        jinds <- add.inds[j] + 1:groups.n[j]
                                        newg <- homophily_cpp$assist$replace_multiple_ties(newg, iinds, jinds, directed, loops)
                                        intergroup.graphs <- rbind(intergroup.graphs, newg)
                                        }
                                    else
                                       intergroup.graphs <- rbind(intergroup.graphs, newg)
                                    }
                                }
                            }
                        }
                    }
                }
            else{
                for (i in 1:(s - 1)){
                    for (j in (i + 1):s){
                        expties <- round(groups.n[i] * groups.n[j] * tie.prob[2])
                        if (expties > 0){
                            send <- add.inds[i] + sample.int(n = groups.n[i], size = expties, replace = TRUE)
                            rec <- add.inds[j] + sample.int(n = groups.n[j], size = expties, replace = TRUE)
                            if (expties == 1)
                                intergroup.graphs <- rbind(intergroup.graphs, c(send, rec))
                            else{
                                newg <- matrix(nrow = expties, ncol = 2)
                                newg[, 1] <- send
                                newg[, 2] <- rec
                                possible.mult.ties <- duplicated(send) * duplicated(rec)
                                if (any(possible.mult.ties)){
                                    iinds <- add.inds[i] + 1:groups.n[i]
                                    jinds <- add.inds[j] + 1:groups.n[j]
                                    newg <- homophily_cpp$assist$replace_multiple_ties(newg, iinds, jinds, directed, loops)
                                    intergroup.graphs <- rbind(intergroup.graphs, newg)
                                    }
                                else
                                    intergroup.graphs <- rbind(intergroup.graphs, newg)
                                }
                            }
                        }
                    }
                }
            g <- as.data.frame(intergroup.graphs)
            g <- graph.data.frame(g, directed = directed, vertices = data.frame(name = 1:n, group = group.assignments))
            }
        else{
            g <- homophily_cpp$games$general(group.assignments, tie.prob, directed, loops)
            g <- graph.adjacency(g, mode = ifelse(directed, "directed", "undirected"), diag = loops)
            }
        }
    V(g)$name <- 1:n
    V(g)$group <- group.assignments
    return(g)
    }
homophily_game <- cmpfun(homophily_game)

group_tie_prob <- function(g, m, mode = "general"){
    # _g_ must be an igraph graph object
    # _m_ must be a vector indicating group membership
    # _mode_ indicates the type of output.
        # "general" returns a vector with two values, one for the density within all groups and the other for between them.
        # "mixture" returns a matrix with densities within and between each group
    if (mode %in% c("general", "mixture") == FALSE)
        stop("_mode_ must be either 'general' or 'mixture'.")
    if (length(m) != vcount(g))
        stop("Length of _m_ does not equal the number of vertices in _g_.")
    g.el <- get.edgelist(g, names = FALSE)
    all.m <- unique(m)
    new.m <- match(m, all.m)
    if (mode == "general")
        out <- homophily_cpp$group_density$general(g.el, new.m, is.directed(g), any(is.loop(g)))
    else{
        out <- homophily_cpp$group_density$mixed(g.el, new.m, is.directed(g), any(is.loop(g)))
        rownames(out) <- as.character(all.m)
        colnames(out) <- as.character(all.m)
        }
    return(out)
    }
group_tie_prob <- cmpfun(group_tie_prob)

###
# In writing the function below, brokerage_roles(), I have found two errors with the comparable function, brokerage(), in library(sna).
# First, brokerage() does not account for partial brokering. See the first paragraph on page 98 of of Gould and Fernandez (1989). brokerage_roles() allows for either implementation.
# Second, for undirected networks, brokerage() uses some double counting.
#     Coordinator, Itinerant broker, and Liaison roles are counted twice.
#     Gatekeeper and Representative seem to be counted correctly.
#     Total brokerage includes the double counted along with adding both Gatekeeper and Representative.
#         (Gatekeeper and Representative are the same and should only count once in the total.)
#     brokerage() is theoretically correct if we think of brokerage as a process of mutual flow in undirected graphs, counting the effect for each brokered vertex.
#     brokerage() is theoretically incorrect if we think of brokerage as a path-based measurement in undirected graphs.
#     Given the emphasis on paths in Gould and Fernandez (1989:95-96), I am inclined to think it's incorrect.
#     (Admittedly, G&F only spoke of directed graphs.)
###

# Returns Coordinator, Itinerant broker, Gatekeeper, Representative, Liaison, and Total.
cppFunction('
    NumericMatrix brokerage_counts(IntegerMatrix twopaths, NumericVector weights, IntegerVector m, bool dir){
        int n = m.size();
        int b;
        int satt, batt, ratt;
        NumericMatrix out(n, 6);
        for (int i = 0; i < n; ++i){
            for (int j = 0; j < 6; ++j){
                out(i, j) = 0.0;
                }
            }
        if (dir){
            for (int i = 0; i < twopaths.nrow(); ++i){
                b = twopaths(i, 1) - 1;
                satt = m[twopaths(i, 0) - 1];
                batt = m[twopaths(i, 1) - 1];
                ratt = m[twopaths(i, 2) - 1];
                if ((satt == batt) & (batt == ratt) & (satt == ratt)){
                    out(b, 0) += weights[i];
                    }
                else if ((satt != batt) & (batt != ratt) & (satt == ratt)){
                    out(b, 1) += weights[i];
                    }
                else if ((satt != batt) & (batt != ratt) & (satt != ratt)){
                    out(b, 4) += weights[i];
                    }
                else if ((satt != batt) & (batt == ratt) & (satt != ratt)){
                    out(b, 2) += weights[i];
                    }
                else{
                    out(b, 3) += weights[i];
                    }
                out(b, 5) += weights[i];
                }
            }
        else{
            for (int i = 0; i < twopaths.nrow(); ++i){
                b = twopaths(i, 1) - 1;
                satt = m[twopaths(i, 0) - 1];
                batt = m[twopaths(i, 1) - 1];
                ratt = m[twopaths(i, 2) - 1];
                if ((satt == batt) & (batt == ratt) & (satt == ratt)){
                    out(b, 0) += weights[i];
                    }
                else if ((satt != batt) & (batt != ratt) & (satt == ratt)){
                    out(b, 1) += weights[i];
                    }
                else if ((satt != batt) & (batt != ratt) & (satt != ratt)){
                    out(b, 4) += weights[i];
                    }
                else{
                    out(b, 2) += weights[i];
                    out(b, 3) += weights[i];
                    }
                out(b, 5) += weights[i];
                }
            }
        return(out);
        }
')

cppFunction('
    IntegerVector group_counts(IntegerVector m){
        int mmax = 0;
        for (int i = 0; i < m.size(); ++i){
            if (m[i] > mmax){
                mmax = m[i];
                }
            }
        IntegerVector out(mmax);
        for (int i = 0; i < mmax; ++i){
            out[i] = 0;
            }
        for (int i = 0; i < m.size(); ++i){
            out[m[i] - 1] += 1;
            }
        return(out);
        }
')

cppFunction('
    NumericMatrix sres2mat(List x, bool partial){
        int nsps = x.size();
        NumericMatrix out(nsps, 4);
        double weight = 1.0 / nsps;
        if (partial){
            for (int i = 0; i < nsps; ++i){
                IntegerVector path = as<IntegerVector>(x[i]);
                out(i, 0) = path[0];
                out(i, 1) = path[1];
                out(i, 2) = path[2];
                out(i, 3) = weight;
                }
            }
        else{
            for (int i = 0; i < nsps; ++i){
                IntegerVector path = as<IntegerVector>(x[i]);
                out(i, 0) = path[0];
                out(i, 1) = path[1];
                out(i, 2) = path[2];
                out(i, 3) = 1.0;
                }
            }
        return(out);
        }
')

brokerage_cpp <- list(brokerage_counts = brokerage_counts, group_counts = group_counts, sp2mat = sres2mat)
rm(brokerage_counts, group_counts, sres2mat)

brokerage_roles <- function(g, m, partial = TRUE, statistics = FALSE, iter = 1000){
    # _g_ must be an igraph graph object
    # _m_ must be a vector indicating group membership
    # _partial_ indicates whether brokerage roles can be shared. (See above comment)
    # If TRUE, _statistics_ performs an attribute permutation test of significance
    # _iter_ are the number of iterations to perform in the permutation tests
    mint <- sort(unique(m))
    mint <- match(m, mint) # m as an integer
    directed <- is.directed(g)
    weightsasone <- rep(1, times = ecount(g)) # Remove edge weights if they exist
    sp <- shortest.paths(g, weights = weightsasone, mode = "out")
    twopathsinds <- which(sp == 2, arr.ind = TRUE)
    twopaths <- lapply(1:nrow(twopathsinds), function(x){
             y <- get.all.shortest.paths(g, twopathsinds[x, 1], twopathsinds[x, 2], mode = "out", weights = weightsasone)$res
             y <- brokerage_cpp$sp2mat(y, partial)
             return(y)
            }
        )
    twopaths <- do.call(rbind, twopaths)
    colnames(twopaths) <- c("sender", "broker", "receiver", "weight")
    if (directed == FALSE){ # Remove duplicate paths
        twopaths <- twopaths[which(twopaths[, "sender"] < twopaths[, "receiver"]), ]
        }

    raw.nli <- brokerage_cpp$brokerage_counts(twopaths[, c("sender", "broker", "receiver")], twopaths[,c("weight")], mint, directed)
    colnames(raw.nli) <- c("Coordinator", "Itinerant Broker", "Gatekeeper", "Representative", "Liaison", "Total")

    raw.gli <- colSums(raw.nli)

    N <- brokerage_cpp$group_counts(mint)
    names(N) <- sort(unique(m))

    if (statistics){
       sims <- replicate(iter, brokerage_cpp$brokerage_counts(twopaths[, c("sender", "broker", "receiver")], twopaths[,c("weight")], sample(mint), directed))
       exp.nli <- apply(sims, c(1,2), mean)[, -6]
       colnames(exp.nli) <- c("Coordinator", "Itinerant Broker", "Gatekeeper", "Representative", "Liaison")
       sd.nli <- apply(sims, c(1,2), sd)[, -6]
       colnames(sd.nli) <- c("Coordinator", "Itinerant Broker", "Gatekeeper", "Representative", "Liaison")
       z.nli <- (raw.nli[, -6] - exp.nli) / sd.nli
       colnames(z.nli) <- c("Coordinator", "Itinerant Broker", "Gatekeeper", "Representative", "Liaison")

       total <- t(apply(sims, 3, colSums))[, -6]
       colnames(total) <- c("Coordinator", "Itinerant Broker", "Gatekeeper", "Representative", "Liaison")
       exp.gli <- colMeans(total)
       sd.gli <- apply(total, 2, sd)
       z.gli <- (raw.gli[-6] - exp.gli) / sd.gli

       return(list(nli = list(raw = raw.nli, exp = exp.nli, sd = sd.nli, z = z.nli), gli = list(raw = raw.gli, exp = exp.gli, sd = sd.gli, z = z.gli), group_counts = N))
       }
    else{
        return(list(nli = raw.nli, gli = raw.gli, group_counts = N))
        }
    }
brokerage_roles <- cmpfun(brokerage_roles)

transitivity_directed <- function(g, method = "weak"){
    # _g_ is an igraph graph object
    # _method_ is the technique of measuring transitivity.
        # "Weak" is the most common understanding of transitivity, so that's the one used here.
        # "Weak census" returns the counts, not the proportion.
        # Once I find a convenient table showing the weights for vacuous transitive triples and see the formula, I'll include strong transitivity.

    if (is.directed(g) == FALSE){
        warning("Entered an undirected graph into transitivity_directed(). Returned transitivity(g).")
        return(transitivity(g))
        }
    acceptable_methods <- c("weak", "weakcensus")
    if (method %in% acceptable_methods == FALSE)
        stop(paste("_method_ must be in c(\"", paste(acceptable_methods, collapse = "\", \""), "\")", sep = ""))
    # See Wasserman and Faust (1994:573)
    weight_mat <- data.frame(row.names = c("003", "012", "102", "021D", "021U", "021C", "111D", "111U", "030T", "030C", "201", "120D", "120U", "120C", "210", "300"))
    weight_mat$b_in <- c(0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 2, 1, 1, 2, 3)
    weight_mat$b_out <- c(0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 2, 1, 2, 3)
    weight_mat$m <- c(0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 2, 1, 1, 1, 2, 3)
    weight_mat$a <- c(0, 1, 0, 2, 2, 2, 1, 1, 3, 3, 0, 2, 2, 2, 1, 0)
    weight_mat$n <- c(3, 2, 2, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0)
    weight_mat$c <- c(0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 6)
    weight_mat$Trans <- c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 2, 1, 3, 6)
    #weight_mat$Vac_trans <- c(6, 6, 6, 6, 6, 5, 5, 5, 5, 3, 4, 4, 4, 3, 2, 0)
    weight_mat$Intrans <- c(0, 0, 0, 0, 0, 1, 1, 1, 0, 3, 2, 0, 0, 2, 1, 0)
    weight_mat$Close_Friends_Disagreeing <- c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 1, 1, 0)

    tc <- triad.census(g)
    if (method == "weakcensus"){
        transcounts <- sum(weight_mat$Trans * tc)
        return(transcounts)
        }
    else if (method == "weak"){
        transcounts <- sum(weight_mat$Trans * tc)
        intranscounts <- sum(weight_mat$Intrans * tc)
        return(transcounts / (transcounts + intranscounts))
        }
    }
transitivity_directed <- cmpfun(transitivity_directed)
# Some additional games that are not in igraph.

# Find a memory-friendly way to do the next two functions without matrices (rescanning the edges.)
cppFunction('
    LogicalMatrix uman_game_prob_cpp(int n, NumericVector dc){
        LogicalMatrix out(n, n);
        int adddyadrandint = 0;
        double adddyadranddoub = 0;
        int addmutrandint = 0;
        double addmutranddoub = 0;
        int dironerandint = 0;
        double dironeranddoub = 0;
        bool adddyad = false;
        bool addmut = false;
        bool dirone = false;
        for (int i = 0; i < (n - 1); ++i){
            for (int j = (i + 1); j < n; ++j){
                adddyadrandint = rand() % 10000000 + 1;
                adddyadranddoub = adddyadrandint / 10000000.0;
                adddyad = adddyadranddoub > dc[2];
                if (adddyad){
                    addmutrandint = rand() % 10000000 + 1;
                    addmutranddoub = addmutrandint / 10000000.0;
                    addmut = addmutranddoub <= (dc[0] / (dc[0] + dc[1]));
                    if (addmut){
                        out(i, j) = true;
                        out(j, i) = true;
                        }
                    else{
                        dironerandint = rand() % 10000000 + 1;
                        dironeranddoub = dironerandint / 10000000.0;
                        dirone = dironeranddoub > .5;
                        if (dirone){
                            out(i, j) = true;
                            out(j, i) = false;
                            }
                        else{
                            out(j, i) = true;
                            out(i, j) = false;
                            }
                        }
                    }
                else{
                    out(i, j) = false;
                    out(j, i) = false;
                    }
                }
            }
        return out;
        }
')

cppFunction('
    IntegerMatrix uman_game_exact_cpp(int n, IntegerVector dc){
        IntegerMatrix out((dc[0] * 2) + dc[1], 2);
        LogicalMatrix usededges(n, n);
        int m = 0;
        int a = 0;
        int counter = 0;
        int r = 0;
        int c = 0;
        int addmutint = 0;
        double addmutdoub = 0.0;
        bool addmut = false;
        for (int i = 0; i < n; ++i){
            for (int j = i; j < n; ++j){
                usededges(i, j) = false;
                usededges(j, i) = false;
                }
            }
        while (counter < out.nrow()){
            r = rand() % n;
            c = rand() % n;
            if ((r != c) & (usededges(r, c) == false)){
                usededges(r, c) = true;
                usededges(c, r) = true;
                if ((m < dc[0]) & (a < dc[1])){
                    addmutint = rand() % 10000000 + 1;
                    addmutdoub = addmutint / 10000000.0;
                    addmut = addmutdoub > (dc[1] / (dc[0] + dc[1]));
                    }
                else if (m == dc[0]){
                    addmut = false;
                    }
                else{
                    addmut = true;
                    }
                if (addmut){
                    out(counter, 0) = r + 1;
                    out(counter, 1) = c + 1;
                    counter += 1;
                    out(counter, 1) = r + 1;
                    out(counter, 0) = c + 1;
                    counter += 1;
                    m += 1;
                    }
                else{
                    out(counter, 0) = r + 1;
                    out(counter, 1) = c + 1;
                    counter += 1;
                    a += 1;
                    }
                }
            }
        return out;
        }
')

cppFunction('
    IntegerMatrix uman_game_exact_sparse_cpp(int n, IntegerVector dc){
        IntegerMatrix out((dc[0] * 2) + dc[1], 2);
        int m = 0;
        int a = 0;
        int counter = 0;
        int r = 0;
        int c = 0;
        int addmutint = 0;
        double addmutdoub = 0.0;
        bool addmut = false;
        while (counter < out.nrow()){
            r = rand() % n;
            c = rand() % n;
            bool used = false;
            int j = 0;
            while ((j < counter) & (used == false)){
                if (((out(j, 0) == (r + 1)) & (out(j, 1) == (c + 1))) | ((out(j, 0) == (c + 1)) & (out(j, 1) == (r + 1)))){
                    used = true;
                    }
                j += 1;
                }
            if ((r != c) & (used == false)){
                if ((m < dc[0]) & (a < dc[1])){
                    addmutint = rand() % 10000000 + 1;
                    addmutdoub = addmutint / 10000000.0;
                    addmut = addmutdoub > (dc[1] / (dc[0] + dc[1]));
                    }
                else if (m == dc[0]){
                    addmut = false;
                    }
                else{
                    addmut = true;
                    }
                if (addmut){
                    out(counter, 0) = r + 1;
                    out(counter, 1) = c + 1;
                    counter += 1;
                    out(counter, 1) = r + 1;
                    out(counter, 0) = c + 1;
                    counter += 1;
                    m += 1;
                    }
                else{
                    out(counter, 0) = r + 1;
                    out(counter, 1) = c + 1;
                    counter += 1;
                    a += 1;
                    }
                }
            }
        return out;
        }
')

.uman_game_cpp <- list(prob = uman_game_prob_cpp, exact = uman_game_exact_cpp, exact_sparse = uman_game_exact_sparse_cpp)
rm(uman_game_prob_cpp, uman_game_exact_cpp, uman_game_exact_sparse_cpp)

uman_game <- function(n, dyad.dist, sparse = FALSE){
    # _n_ is the number of vertices in the desired graph
    # _dyad.dist_ is a vector of length three. The values correspond to Mutual, Asymmetric, and Null dyads
      # If the values are all less than or equal to 1 (and greater than 0), they're probabilities
      # If any of the values are greater than 1, they're exact edge counts
    # _sparse_ is good for large networks with low densities. It's poor for dense and small networks.
    n.dyads <- n * (n - 1) / 2
    if (any(dyad.dist > 1) & (sum(round(dyad.dist)) != n.dyads))
        stop("Number of dyads is not possible given the number of vertices. (Multiple edges are not supported.)")
    if (length(dyad.dist) != 3)
        stop("_dyad.dist_ corresponds to the dyad counts or probabilities. Its length must equal three, corresponding to mutual, asymmetric, and null dyads (in that order).")

    if (all(dyad.dist <= 1)){
        if (sparse){
            dyad.dist.round <- round(dyad.dist * n.dyads)
            dyad.dist.diff <- sum(dyad.dist.round) - n.dyads
            # If there's a rounding difference, adjust the nulls.
            dyad.dist.round[3] <- dyad.dist.round[3] - dyad.dist.diff
            g <- as.data.frame(.uman_game_cpp$exact_sparse(n, dyad.dist.round))
            g <- graph.data.frame(g, directed = TRUE, vertices = data.frame(name = 1:n))
            }
        else{
            dyad.dist <- round(dyad.dist, digits = 6)
            g <- .uman_game_cpp$prob(n, dyad.dist)
            g <- graph.adjacency(g, mode = "directed", diag = FALSE)
            }
        }
    else{
        dyad.dist <- round(dyad.dist)
        if (sparse)
            g <- .uman_game_cpp$exact_sparse(n, dyad.dist)
        else
            g <- .uman_game_cpp$exact(n, dyad.dist)
        g <- as.data.frame(g)
        g <- graph.data.frame(g, directed = TRUE, vertices = data.frame(name = 1:n))
        }
    return(g)
    }
uman_game <- cmpfun(uman_game)

# Rewires the dyads. It creates a random graph with the same:
#     Indegree and outdegree for each vertex (and density)
#     The level of reciprocity for each vertex (and dyad census)
# It will not produce multiple edges or loops
# I would not recommend this function for very large networks, as the matrix will get big.
# Later, if needed, I can create a sparse version that stores ties as an adjacency list rather than a matrix.

cppFunction('
    IntegerMatrix dyad_rewire_cpp(IntegerMatrix el, int nswaps){
        int ne = el.nrow();
        IntegerMatrix out(ne, 2);
        int i, j, asymcount, nd, na, d1, d2, s1, r1, s2, r2;
        int nv = 0;
        int counter = 0;
        int nmut = 0;
        int mutcount = 0;
        bool diffind, nonloop, mutswap, nonexist, nonrecipexist;

        for (i = 0; i < ne; ++i){
            if (el(i, 0) > nv){
                nv = el(i, 0);
                }
            if (el(i, 1) > nv){
                nv = el(i, 1);
                }
            }
        LogicalMatrix mat(nv, nv);
        for (i = 0; i < (nv - 1); ++i){
            mat(i, i) = false;
            for (j = (i + 1); j < nv; ++j){
                mat(i, j) = false;
                mat(j, i) = false;
                }
            }
        for (i = 0; i < ne; ++i){
            mat((el(i, 0) - 1), (el(i, 1) - 1)) = true;
            if (mat((el(i, 0) - 1), (el(i, 1) - 1)) & mat((el(i, 1) - 1), (el(i, 0) - 1))){
                nmut += 1;
                }
            }
        nd = ne - nmut;
        na = nd - nmut;
        IntegerMatrix dl(nd, 2);
        asymcount = nmut;
        for (i = 0; i < (nv - 1); ++i){
            for (j = (i + 1); j < nv; ++j){
                if (mat(i, j) & mat(j, i)){
                    dl(mutcount, 0) = i + 1;
                    dl(mutcount, 1) = j + 1;
                    mutcount += 1;
                    }
                else if (mat(i, j) & (mat(j, i) == false)){
                    dl(asymcount, 0) = i + 1;
                    dl(asymcount, 1) = j + 1;
                    asymcount += 1;
                    }
                else if (mat(j, i) & (mat(i, j) == false)){
                    dl(asymcount, 0) = j + 1;
                    dl(asymcount, 1) = i + 1;
                    asymcount += 1;
                    }
                }
            }

        while (counter < nswaps){
            d1 = rand() % nd;
            mutswap = d1 < nmut;
            if (mutswap){
                d2 = rand() % nmut;
                }
            else{
                d2 = rand() % na;
                d2 += nmut;
                }
            s1 = dl(d1, 0);
            r1 = dl(d1, 1);
            s2 = dl(d2, 0);
            r2 = dl(d2, 1);
            diffind = (d1 != d2);
            nonloop = ((s1 != r2) & (s2 != r1));
            nonexist = ((mat(s1 - 1, r2 - 1) == false) & (mat(s2 - 1, r1 - 1) == false));
            nonrecipexist = ((mat(r2 - 1, s1 - 1) == false) & (mat(r1 - 1, s2 - 1) == false));
            if (diffind & nonloop & nonexist & nonrecipexist){
                if (mutswap){
                    dl(d1, 1) = r2;
                    dl(d2, 1) = r1;
                    mat(s1 - 1, r2 - 1) = true;
                    mat(s1 - 1, r1 - 1) = false;
                    mat(r2 - 1, s1 - 1) = true;
                    mat(r1 - 1, s1 - 1) = false;

                    mat(s2 - 1, r1 - 1) = true;
                    mat(s2 - 1, r2 - 1) = false;
                    mat(r1 - 1, s2 - 1) = true;
                    mat(r2 - 1, s2 - 1) = false;
                    counter += 1;
                    }
                else{
                    dl(d1, 1) = r2;
                    dl(d2, 1) = r1;
                    mat(s1 - 1, r2 - 1) = true;
                    mat(s1 - 1, r1 - 1) = false;
                    mat(s2 - 1, r1 - 1) = true;
                    mat(s2 - 1, r2 - 1) = false;
                    counter += 1;
                    }
                }
            }
        counter = 0;
        for (i = 0; i < nmut; ++i){
            out(counter, 0) = dl(i, 0);
            out(counter, 1) = dl(i, 1);
            counter += 1;
            out(counter, 0) = dl(i, 1);
            out(counter, 1) = dl(i, 0);
            counter += 1;
            }
        for (i = nmut; i < nd; ++i){
            out(counter, 0) = dl(i, 0);
            out(counter, 1) = dl(i, 1);
            counter += 1;
            }
        return out;
        }
')

.dyad_rewire_cpp <- dyad_rewire_cpp
rm(dyad_rewire_cpp)

rewire_dyads <- function(g, niter = 100){
    # _g_ is a directed igraph object
    # _niter_ is the number of dyad rewires to perform.

    if (is.directed(g) == FALSE)
        stop("_g_ should be a directed igraph object.")
    if (is.simple(g) == FALSE)
        stop("Graph _g_ must be simple. Multiple edges and loops are not currently supported.")

    gel <- get.edgelist(g, names = FALSE)
    gel <- .dyad_rewire_cpp(gel, niter)
    gel <- as.data.frame(gel)
    # Perhaps consider adding in the vertex and edge attributes here.
    newg <- graph.data.frame(gel, directed = TRUE, vertices = data.frame(name = 1:vcount(g)))
    return(newg)
    }
rewire_dyads <- cmpfun(rewire_dyads)

# See David Krackhardt. 1994. "Graph Theoretical Dimensions of Informal Organizations" in Kathleen Carley and Michael Prietula (eds.) _Computational Organizational Theory._ Hillsdale, NJ: Lawrence Earlbaum Associates, Inc, p. 89 - 111.

cppFunction('
    int lubness_con_cpp(LogicalMatrix r){
        int i, j, k, l, n, nub, viol;
        char lub, notlub;
        viol = 0;
        n = r.nrow();
        IntegerVector ub(n);

        if (n > 2){
            for (i = 0; i < (n - 1); i++){
                for (j = i + 1; j < n; j++){
                    nub = 0;
                    for (k = 0; k < n; k++){
                        if (r(k, i) && r(k, j)){
                            ub[nub] = k;
                            nub += 1;
                            }
                        }
                    lub = 0;
                    for (k = 0; (k < nub) && (!lub); k++){
                        notlub = 0;
                        for (l = 0; (l < nub) && (!notlub); l++){
                            if (!r(ub[k], ub[l])){
                                notlub += 1;
                                }
                            }
                        if (!notlub){
                            lub += 1;
                            }
                        }
                    if (!lub){
                        viol += 1;
                        }
                    }
                }
            }
        return viol;
        }
')

cppFunction('
    IntegerVector V_MaxV_Effic(int Nn, int Ne){
        IntegerVector out(2);
        out[0] = Ne - (Nn - 1);
        out[1] = (Nn * (Nn - 1)) - (Nn - 1);
        return out;
        }
')

cppFunction('
    IntegerVector rep1(int N){
        IntegerVector out(N);
        for (int i = 0; i < N; i++)
            out[i] = 1;
        return out;
        }
')

cppFunction('
    IntegerVector repx(int x, int N){
        IntegerVector out(N);
        for (int i = 0; i < N; i++)
            out[i] = x;
        return out;
        }
')

cppFunction('
    int sumTRUE(LogicalVector x){
        int out = 0;
        int n = x.size();

        for (int i = 0; i < n; i++)
            if (x[i])
                out += 1;
        return out;
        }
')

cppFunction('
    int intsum(IntegerVector x){
        int out = 0;
        int n = x.size();

        for (int i = 0; i < n; i++)
            out += x[i];
        return out;
        }
')

.krackhardt_1994 <- list(lubness_con = lubness_con_cpp, v_maxv_effic = V_MaxV_Effic, utils = list(rep1 = rep1, repx = repx, sumTRUE = sumTRUE, intsum = intsum))

rm(lubness_con_cpp, V_MaxV_Effic, rep1, repx, sumTRUE, intsum)

connectedness <- function(g){
    if (is.directed(g) == FALSE)
        stop("Graph _g_ must be directed.")
    if (any(is.loop(g)))
        stop("Graph _g_ must not have loops.")
    N <- vcount(g)
    allpairs <- N * (N - 1) / 2
    eweight <- .krackhardt_1994$utils$rep1(ecount(g))
    discons <- sapply(1:(N - 1), function(x){
            targets <- (x + 1):N
            temppls <- shortest.paths(g, v = x, to = targets, mode = "all", weights = eweight)
            discon <- is.infinite(temppls)
            discon <- .krackhardt_1994$utils$sumTRUE(discon)
            return(discon)
            }
        )
    return(1 - (.krackhardt_1994$utils$intsum(discons) / allpairs))
    }
connectedness <- cmpfun(connectedness)

hierarchy_krackhardt <- function(g){
    if (is.directed(g) == FALSE)
        stop("Graph _g_ must be directed.")
    if (any(is.loop(g)))
        stop("Graph _g_ must not have loops.")
    N <- vcount(g)
    eweight <- .krackhardt_1994$utils$rep1(ecount(g))
    conned <- lapply(1:N, function(x){
            targets <- c(1:N)[-x]
            temppls <- shortest.paths(g, v = x, to = targets, mode = "out", weights = eweight)
            reachable <- targets[which(is.finite(temppls))]
            reachableedges <- cbind(rep(x, times = length(reachable)), reachable)
            return(reachableedges)
            }
        )
    conned <- do.call(rbind, conned)
    conned <- graph.edgelist(as.matrix(conned), directed = TRUE)
    conneddc <- unlist(dyad.census(conned))
    return(1 - (conneddc["mut"] / (conneddc["mut"] + conneddc["asym"])))
    }
hierarchy_krackhardt <- cmpfun(hierarchy_krackhardt)

# The code for efficiency runs ~3x slower than the library(sna) code on the small EIES data.
efficiency <- function(g){
    if (is.directed(g) == FALSE)
        stop("Graph _g_ must be directed.")
    if (is.simple(g) == FALSE)
        stop("Graph _g_ must be a simple graph.")
    gclust <- clusters(g, mode = "weak")
    excesscount <- lapply(1:gclust$no, function(x){
            subg <- induced.subgraph(g, vids = which(gclust$membership == x))
            Ne <- ecount(subg)
            Nn <- vcount(subg)
            return(.krackhardt_1994$v_maxv_effic(Nn, Ne))
            }
        )
    excesscount <- do.call(rbind, excesscount)
    excesscount <- colSums(excesscount) # Wrote an Rcpp colSums equivalent. R's native colSums works faster.
    return(1 - (excesscount[1] / excesscount[2]))
    }
efficiency <- cmpfun(efficiency)

# The LUBness program was initially written by Carter T. Butts for library(sna). I have adapted the code somewhat for igraph and Rcpp.

lubness <- function(g){
    if (is.directed(g) == FALSE)
        stop("Graph _g_ must be directed.")
    if (any(is.loop(g)))
        stop("Graph _g_ must not have loops.")
    cd <- clusters(g, mode = "weak")
    nolub <- 0
    maxnolub <- 0
    for (i in 1:cd$no){
        vi <- which(cd$membership == i)
        subg <- induced.subgraph(g, vi)
        ne <- ecount(subg)
        E(subg)$weight <- rep(1, times = ne) # Use Rcpp here, rep1
        nv <- vcount(subg)
        r <- lapply(1:nv, function(x){
                sp <- shortest.paths(subg, v = x, mode = "out")
                sp <- is.finite(sp)
                return(sp)
                }
            )
        r <- do.call(rbind, r)

        if (cd$csize[i] > 2){
            viol <- .krackhardt_1994$lubness_con(r)
            nolub <- nolub + viol
            maxnolub <- maxnolub + (cd$csize[i] - 1) * (cd$csize[i] - 2) / 2 # Use Rcpp here
            }
        }
    return(1 - nolub / maxnolub)
    }
lubness <- cmpfun(lubness)

# My C++ code generally performs poorer than my R code
vazquez_game<-function(N, u){
  # _N_ is the number of nodes
  # _u_ is the triadic closure parameter. Needs to be between 0 and 1.
  # Alexei Vazquez.  2003.  "Growing network with local rules: Preferential attachment, clustering hierarchy, and degree correlations."  Physical Review E 67, 056104
  add.node.edge<-function()
    return(add.node = c(rbinom(1, 1, 1 - u), close.tri = rbinom(1, 1, u)))
  vneigh <- function(v, el){
    receivers <- el[which(el[ , 1] == v), 2]
    senders <- el[which(el[ , 2] == v), 1]
    return(unique(c(receivers, senders)))
    }
  g <- matrix(NA, 0, 2)
  potentials <- matrix(NA, 0, 2)
  ncount <- 1
  while(ncount < N){
    a <- add.node.edge()
    if (a[1] == 1){
      source.v <- ncount + 1
      target.v <- ifelse(ncount > 1, sample.int(ncount, 1), 1)
      target.neighbors <- vneigh(target.v, g)
      g <- rbind(g, c(source.v, target.v))
      n.neighbors <- length(target.neighbors)
      if (n.neighbors > 0){
        new.potentials <- matrix(NA, n.neighbors, 2)
        new.potentials[ , 1] <- source.v
        new.potentials[ , 2] <- target.neighbors
        potentials <- rbind(potentials, new.potentials)
        rm(new.potentials)
        }
      rm(source.v, target.v, target.neighbors, n.neighbors)
      ncount <- ncount + 1
      }
    if ((a[2] == 1) & (ncount > 2) & (length(potentials) >= 2)){
      if (length(potentials) == 2){
        g <- rbind(g, c(potentials[1], potentials[2]))
        potentials <- matrix(NA, 0, 2)
        }
      else{
        potentials.ind <- sample.int(nrow(potentials), 1)
        adders <- potentials[potentials.ind, ]
        g <- rbind(g, c(adders[1], adders[2]))
        potentials <- potentials[-potentials.ind, ]
        }
      }
    }
  g <- graph.edgelist(g, directed = FALSE)
  return(g)
  }
vazquez_game<-cmpfun(vazquez_game)