API tools faq
paste
BenjaminLind

AdditionalFunctionsIgraph.R

BenjaminLind
Aug 8th, 2014
531
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 16.67 KB | None | 0 0
raw download clone embed print report
  1. #####
  2. ### Except where otherwise indicated, the code below is written by Benjamin Lind under CC: BY-NC-SA 3.0
  3. ### If you use or develop this code, please provide proper attribution.
  4. ### You may contact me at lind.benjamin//at//gmail//d0t//com
  5. ### You may also find me @benjaminlind on Twitter
  6. #####
  7.  
  8. ### To load these functions into R, enter:
  9. ###     source("http://pastebin.com/raw.php?i=SErkr5sg")
  10.  
  11. ### Functions included in this file:
  12. ###     Krackhardt's E-I index
  13. ###     Weak transitivity as measured for a directed network
  14. ###     UMAN random graph generators. They can make random graphs conditioned on the dyad census.
  15. ###     A dyad rewired graph. Like igraph::rewire(), but controls for reciprocity.
  16.  
  17. library(Rcpp)
  18. library(igraph)
  19. library(compiler)
  20.  
  21. cppFunction('
  22.   double ei_index_cpp(IntegerMatrix el, IntegerVector attrib){
  23.       double incount = 0;
  24.       int n = el.nrow();
  25.       for (int i = 0; i < n; ++i){
  26.           int attrib1 = attrib[el(i, 0) - 1];
  27.           int attrib2 = attrib[el(i, 1) - 1];
  28.           if (attrib1 == attrib2){
  29.               incount += 1;
  30.               }
  31.           }
  32.       double excount = n - incount;
  33.       double out = (excount - incount) / n;
  34.       return out;
  35.       }
  36. ')
  37.  
  38. ei_index <- function(g, value){
  39.     # _g_ is an igraph 'graph' object
  40.     # _value_ is a categorical vertex attribute
  41.     if (any(is.na(value)))
  42.         stop("Missing values are not allowed for _value_.")
  43.     if (vcount(g) != length(value))
  44.         stop("The length of _value_ does not equal the number of vertices in _g_.")
  45.     g.el <- get.edgelist(g, names = FALSE)
  46.     allvalues <- unique(value)
  47.     value.int <- match(value, allvalues)
  48.     return(ei_index_cpp(g.el, value.int))
  49.     }
  50. ei_index <- cmpfun(ei_index)
  51.  
  52. transitivity_directed <- function(g, method = "weak"){
  53.     # _g_ is an igraph graph object
  54.     # _method_ is the technique of measuring transitivity.
  55.         # "Weak" is the most common understanding of transitivity, so that's the one used here.
  56.         # "Weak census" returns the counts, not the proportion.
  57.         # Once I find a convenient table showing the weights for vacuous transitive triples and see the formula, I'll include strong transitivity.
  58.  
  59.     if (is.directed(g) == FALSE){
  60.         warning("Entered an undirected graph into transitivity_directed(). Returned transitivity(g).")
  61.         return(transitivity(g))
  62.         }    
  63.     acceptable_methods <- c("weak", "weakcensus")
  64.     if (method %in% acceptable_methods == FALSE)
  65.         stop(paste("_method_ must be in c(\"", paste(acceptable_methods, collapse = "\", \""), "\")", sep = ""))
  66.     # See Wasserman and Faust (1994:573)
  67.     weight_mat <- data.frame(row.names = c("003", "012", "102", "021D", "021U", "021C", "111D", "111U", "030T", "030C", "201", "120D", "120U", "120C", "210", "300"))
  68.     weight_mat$b_in <- c(0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 2, 1, 1, 2, 3)
  69.     weight_mat$b_out <- c(0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 2, 1, 2, 3)
  70.     weight_mat$m <- c(0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 2, 1, 1, 1, 2, 3)
  71.     weight_mat$a <- c(0, 1, 0, 2, 2, 2, 1, 1, 3, 3, 0, 2, 2, 2, 1, 0)
  72.     weight_mat$n <- c(3, 2, 2, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0)
  73.     weight_mat$c <- c(0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 6)
  74.     weight_mat$Trans <- c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 2, 1, 3, 6)
  75.     #weight_mat$Vac_trans <- c(6, 6, 6, 6, 6, 5, 5, 5, 5, 3, 4, 4, 4, 3, 2, 0)
  76.     weight_mat$Intrans <- c(0, 0, 0, 0, 0, 1, 1, 1, 0, 3, 2, 0, 0, 2, 1, 0)
  77.     weight_mat$Close_Friends_Disagreeing <- c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 1, 1, 0)
  78.  
  79.     tc <- triad.census(g)
  80.     if (method == "weakcensus"){
  81.         transcounts <- sum(weight_mat$Trans * tc)
  82.         return(transcounts)
  83.         }
  84.     else if (method == "weak"){
  85.         transcounts <- sum(weight_mat$Trans * tc)
  86.         intranscounts <- sum(weight_mat$Intrans * tc)
  87.         return(transcounts / (transcounts + intranscounts))
  88.         }
  89.     }
  90. transitivity_directed <- cmpfun(transitivity_directed)
  91.  
  92. cppFunction('    
  93.   LogicalMatrix uman_game_prob_cpp(int n, NumericVector dc){
  94.       LogicalMatrix out(n, n);
  95.       int adddyadrandint = 0;
  96.       double adddyadranddoub = 0;
  97.       int addmutrandint = 0;
  98.       double addmutranddoub = 0;
  99.       int dironerandint = 0;
  100.       double dironeranddoub = 0;
  101.       bool adddyad = false;
  102.       bool addmut = false;
  103.       bool dirone = false;
  104.       for (int i = 0; i < (n - 1); ++i){
  105.           for (int j = (i + 1); j < n; ++j){
  106.               adddyadrandint = rand() % 10000000 + 1;
  107.               adddyadranddoub = adddyadrandint / 10000000.0;
  108.               adddyad = adddyadranddoub > dc[2];
  109.               if (adddyad){
  110.                   addmutrandint = rand() % 10000000 + 1;
  111.                   addmutranddoub = addmutrandint / 10000000.0;
  112.                   addmut = addmutranddoub <= (dc[0] / (dc[0] + dc[1]));
  113.                   if (addmut){
  114.                       out(i, j) = true;
  115.                       out(j, i) = true;
  116.                       }
  117.                   else{
  118.                       dironerandint = rand() % 10000000 + 1;
  119.                       dironeranddoub = dironerandint / 10000000.0;
  120.                       dirone = dironeranddoub > .5;
  121.                       if (dirone){
  122.                           out(i, j) = true;
  123.                           out(j, i) = false;
  124.                           }
  125.                       else{
  126.                           out(j, i) = true;
  127.                           out(i, j) = false;
  128.                           }
  129.                       }
  130.                   }
  131.               else{
  132.                   out(i, j) = false;
  133.                   out(j, i) = false;                    
  134.                   }
  135.               }
  136.           }
  137.       return out;
  138.       }
  139. ')
  140.  
  141. cppFunction('    
  142.   IntegerMatrix uman_game_exact_cpp(int n, IntegerVector dc){
  143.       IntegerMatrix out((dc[0] * 2) + dc[1], 2);
  144.       LogicalMatrix usededges(n, n);
  145.       int m = 0;
  146.       int a = 0;
  147.       int counter = 0;
  148.       int r = 0;
  149.       int c = 0;
  150.       int addmutint = 0;
  151.       double addmutdoub = 0.0;
  152.       bool addmut = false;
  153.       for (int i = 0; i < n; ++i){
  154.           for (int j = i; j < n; ++j){
  155.               usededges(i, j) = false;
  156.               usededges(j, i) = false;
  157.               }
  158.           }
  159.       while (counter < out.nrow()){
  160.           r = rand() % n;
  161.           c = rand() % n;
  162.           if ((r != c) & (usededges(r, c) == false)){
  163.               usededges(r, c) = true;
  164.               usededges(c, r) = true;
  165.               if ((m < dc[0]) & (a < dc[1])){
  166.                   addmutint = rand() % 10000000 + 1;
  167.                   addmutdoub = addmutint / 10000000.0;
  168.                   addmut = addmutdoub > (dc[1] / (dc[0] + dc[1]));
  169.                   }
  170.               else if (m == dc[0]){
  171.                   addmut = false;
  172.                   }
  173.               else{
  174.                   addmut = true;
  175.                   }
  176.               if (addmut){
  177.                   out(counter, 0) = r + 1;
  178.                   out(counter, 1) = c + 1;
  179.                   counter += 1;
  180.                   out(counter, 1) = r + 1;
  181.                   out(counter, 0) = c + 1;
  182.                   counter += 1;
  183.                   m += 1;
  184.                   }
  185.               else{
  186.                   out(counter, 0) = r + 1;
  187.                   out(counter, 1) = c + 1;
  188.                   counter += 1;                    
  189.                   a += 1;
  190.                   }
  191.               }
  192.           }
  193.       return out;
  194.       }
  195. ')
  196.  
  197. cppFunction('    
  198.   IntegerMatrix uman_game_exact_sparse_cpp(int n, IntegerVector dc){
  199.       IntegerMatrix out((dc[0] * 2) + dc[1], 2);
  200.       int m = 0;
  201.       int a = 0;
  202.       int counter = 0;
  203.       int r = 0;
  204.       int c = 0;
  205.       int addmutint = 0;
  206.       double addmutdoub = 0.0;
  207.       bool addmut = false;
  208.       while (counter < out.nrow()){
  209.           r = rand() % n;
  210.           c = rand() % n;
  211.           bool used = false;
  212.           int j = 0;
  213.           while ((j < counter) & (used == false)){
  214.               if (((out(j, 0) == (r + 1)) & (out(j, 1) == (c + 1))) | ((out(j, 0) == (c + 1)) & (out(j, 1) == (r + 1)))){
  215.                   used = true;
  216.                   }
  217.               j += 1;
  218.               }
  219.           if ((r != c) & (used == false)){
  220.               if ((m < dc[0]) & (a < dc[1])){
  221.                   addmutint = rand() % 10000000 + 1;
  222.                   addmutdoub = addmutint / 10000000.0;
  223.                   addmut = addmutdoub > (dc[1] / (dc[0] + dc[1]));
  224.                   }
  225.               else if (m == dc[0]){
  226.                   addmut = false;
  227.                   }
  228.               else{
  229.                   addmut = true;
  230.                   }
  231.               if (addmut){
  232.                   out(counter, 0) = r + 1;
  233.                   out(counter, 1) = c + 1;
  234.                   counter += 1;
  235.                   out(counter, 1) = r + 1;
  236.                   out(counter, 0) = c + 1;
  237.                   counter += 1;
  238.                   m += 1;
  239.                   }
  240.               else{
  241.                   out(counter, 0) = r + 1;
  242.                   out(counter, 1) = c + 1;
  243.                   counter += 1;                    
  244.                   a += 1;
  245.                   }
  246.               }
  247.           }
  248.       return out;
  249.       }
  250. ')
  251.  
  252. .uman_game_cpp <- list(prob = uman_game_prob_cpp, exact = uman_game_exact_cpp, exact_sparse = uman_game_exact_sparse_cpp)
  253. rm(uman_game_prob_cpp, uman_game_exact_cpp, uman_game_exact_sparse_cpp)
  254.  
  255. uman_game <- function(n, dyad.dist, sparse = FALSE){
  256.     # _n_ is the number of vertices in the desired graph
  257.     # _dyad.dist_ is a vector of length three. The values correspond to Mutual, Asymmetric, and Null dyads
  258.       # If the values are all less than or equal to 1 (and greater than 0), they're probabilities
  259.       # If any of the values are greater than 1, they're exact edge counts
  260.     # _sparse_ is good for large networks with low densities. It's poor for dense and small networks.
  261.     n.dyads <- n * (n - 1) / 2
  262.     if (any(dyad.dist > 1) & (sum(round(dyad.dist)) != n.dyads))
  263.         stop("Number of dyads is not possible given the number of vertices. (Multiple edges are not supported.)")
  264.     if (length(dyad.dist) != 3)
  265.         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).")
  266.  
  267.     if (all(dyad.dist <= 1)){
  268.         if (sparse){
  269.             dyad.dist.round <- round(dyad.dist * n.dyads)
  270.             dyad.dist.diff <- sum(dyad.dist.round) - n.dyads
  271.             # If there's a rounding difference, adjust the nulls.
  272.             dyad.dist.round[3] <- dyad.dist.round[3] - dyad.dist.diff
  273.             g <- as.data.frame(.uman_game_cpp$exact_sparse(n, dyad.dist.round))
  274.             g <- graph.data.frame(g, directed = TRUE, vertices = data.frame(name = 1:n))
  275.             }
  276.         else{
  277.             dyad.dist <- round(dyad.dist, digits = 6)
  278.             g <- .uman_game_cpp$prob(n, dyad.dist)
  279.             g <- graph.adjacency(g, mode = "directed", diag = FALSE)
  280.             }
  281.         }
  282.     else{
  283.         dyad.dist <- round(dyad.dist)
  284.         if (sparse)
  285.             g <- .uman_game_cpp$exact_sparse(n, dyad.dist)
  286.         else
  287.             g <- .uman_game_cpp$exact(n, dyad.dist)
  288.         g <- as.data.frame(g)
  289.         g <- graph.data.frame(g, directed = TRUE, vertices = data.frame(name = 1:n))
  290.         }
  291.     return(g)
  292.     }
  293. uman_game <- cmpfun(uman_game)
  294.  
  295. # Rewires the dyads. It creates a random graph with the same:
  296. #     Indegree and outdegree for each vertex (and density)
  297. #     The level of reciprocity for each vertex (and dyad census)
  298. # It will not produce multiple edges or loops
  299. # I would not recommend this function for very large networks, as the matrix will get big.
  300. # Later, if needed, I can create a sparse version that stores ties as an adjacency list rather than a matrix.
  301.  
  302. cppFunction('
  303.   IntegerMatrix dyad_rewire_cpp(IntegerMatrix el, int nswaps){
  304.       int ne = el.nrow();
  305.       IntegerMatrix out(ne, 2);
  306.       int i, j, asymcount, nd, na, d1, d2, s1, r1, s2, r2;
  307.       int nv = 0;
  308.       int counter = 0;
  309.       int nmut = 0;
  310.       int mutcount = 0;
  311.       bool diffind, nonloop, mutswap, nonexist, nonrecipexist;
  312.              
  313.       for (i = 0; i < ne; ++i){
  314.           if (el(i, 0) > nv){
  315.               nv = el(i, 0);
  316.               }
  317.           if (el(i, 1) > nv){
  318.               nv = el(i, 1);
  319.               }
  320.           }
  321.       LogicalMatrix mat(nv, nv);
  322.       for (i = 0; i < (nv - 1); ++i){
  323.           mat(i, i) = false;
  324.           for (j = (i + 1); j < nv; ++j){
  325.               mat(i, j) = false;
  326.               mat(j, i) = false;
  327.               }
  328.           }
  329.       for (i = 0; i < ne; ++i){            
  330.           mat((el(i, 0) - 1), (el(i, 1) - 1)) = true;
  331.           if (mat((el(i, 0) - 1), (el(i, 1) - 1)) & mat((el(i, 1) - 1), (el(i, 0) - 1))){
  332.               nmut += 1;
  333.               }
  334.           }
  335.       nd = ne - nmut;
  336.       na = nd - nmut;
  337.       IntegerMatrix dl(nd, 2);
  338.       asymcount = nmut;
  339.       for (i = 0; i < (nv - 1); ++i){
  340.           for (j = (i + 1); j < nv; ++j){
  341.               if (mat(i, j) & mat(j, i)){
  342.                   dl(mutcount, 0) = i + 1;
  343.                   dl(mutcount, 1) = j + 1;
  344.                   mutcount += 1;
  345.                   }
  346.               else if (mat(i, j) & (mat(j, i) == false)){
  347.                   dl(asymcount, 0) = i + 1;
  348.                   dl(asymcount, 1) = j + 1;
  349.                   asymcount += 1;
  350.                   }
  351.               else if (mat(j, i) & (mat(i, j) == false)){
  352.                   dl(asymcount, 0) = j + 1;
  353.                   dl(asymcount, 1) = i + 1;
  354.                   asymcount += 1;
  355.                   }
  356.               }
  357.           }      
  358.      
  359.       while (counter < nswaps){
  360.           d1 = rand() % nd;
  361.           mutswap = d1 < nmut;
  362.           if (mutswap){
  363.               d2 = rand() % nmut;
  364.               }
  365.           else{
  366.               d2 = rand() % na;
  367.               d2 += nmut;
  368.               }
  369.           s1 = dl(d1, 0);
  370.           r1 = dl(d1, 1);
  371.           s2 = dl(d2, 0);
  372.           r2 = dl(d2, 1);
  373.           diffind = (d1 != d2);
  374.           nonloop = ((s1 != r2) & (s2 != r1));
  375.           nonexist = ((mat(s1 - 1, r2 - 1) == false) & (mat(s2 - 1, r1 - 1) == false));
  376.           nonrecipexist = ((mat(r2 - 1, s1 - 1) == false) & (mat(r1 - 1, s2 - 1) == false));
  377.           if (diffind & nonloop & nonexist & nonrecipexist){            
  378.               if (mutswap){
  379.                   dl(d1, 1) = r2;
  380.                   dl(d2, 1) = r1;
  381.                   mat(s1 - 1, r2 - 1) = true;
  382.                   mat(s1 - 1, r1 - 1) = false;
  383.                   mat(r2 - 1, s1 - 1) = true;
  384.                   mat(r1 - 1, s1 - 1) = false;
  385.                  
  386.                   mat(s2 - 1, r1 - 1) = true;
  387.                   mat(s2 - 1, r2 - 1) = false;
  388.                   mat(r1 - 1, s2 - 1) = true;
  389.                   mat(r2 - 1, s2 - 1) = false;
  390.                   counter += 1;
  391.                   }
  392.               else{
  393.                   dl(d1, 1) = r2;
  394.                   dl(d2, 1) = r1;
  395.                   mat(s1 - 1, r2 - 1) = true;
  396.                   mat(s1 - 1, r1 - 1) = false;
  397.                   mat(s2 - 1, r1 - 1) = true;
  398.                   mat(s2 - 1, r2 - 1) = false;
  399.                   counter += 1;
  400.                   }
  401.               }
  402.           }
  403.       counter = 0;
  404.       for (i = 0; i < nmut; ++i){
  405.           out(counter, 0) = dl(i, 0);
  406.           out(counter, 1) = dl(i, 1);
  407.           counter += 1;
  408.           out(counter, 0) = dl(i, 1);
  409.           out(counter, 1) = dl(i, 0);
  410.           counter += 1;
  411.           }
  412.       for (i = nmut; i < nd; ++i){
  413.           out(counter, 0) = dl(i, 0);
  414.           out(counter, 1) = dl(i, 1);
  415.           counter += 1;
  416.           }
  417.       return out;
  418.       }
  419. ')
  420.  
  421. .dyad_rewire_cpp <- dyad_rewire_cpp
  422. rm(dyad_rewire_cpp)
  423.  
  424. rewire_dyads <- function(g, niter = 100){
  425.     # _g_ is a directed igraph object
  426.     # _niter_ is the number of dyad rewires to perform.
  427.  
  428.     if (is.directed(g) == FALSE)
  429.         stop("_g_ should be a directed igraph object.")
  430.     if (is.simple(g) == FALSE)
  431.         stop("Graph _g_ must be simple. Multiple edges and loops are not currently supported.")
  432.  
  433.     gel <- get.edgelist(g, names = FALSE)
  434.     gel <- .dyad_rewire_cpp(gel, niter)
  435.     gel <- as.data.frame(gel)
  436.     # Perhaps consider adding in the vertex and edge attributes here.
  437.     newg <- graph.data.frame(gel, directed = TRUE, vertices = data.frame(name = 1:vcount(g)))
  438.     return(newg)
  439.     }
  440. rewire_dyads <- cmpfun(rewire_dyads)
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!