PAFit: Generative Mechanism Estimation in Temporal Complex Networks

Documented in plot.PAFit_result

# function to plot estimation results  2015-3-12 Thong Pham
plot.PAFit_result <-
function(x                       ,
         net_stat       = NULL   ,
         true_f         = NULL   , plot             = "A"              , plot_bin     = TRUE  ,
         line           = FALSE  , confidence       = TRUE             , high_deg_A   = 1     , 
         high_deg_f     = 5      ,   
         shade_point    = 0.5    , col_point        = "grey25"         , pch          = 16    , 
         shade_interval = 0.5    , col_interval     = "lightsteelblue" ,
         label_x        = NULL   , label_y          = NULL             ,
         max_A          = NULL   , min_A            = NULL             , f_min        = NULL  , 
         f_max          = NULL   , plot_true_degree = FALSE            , 
         ...) {
  if (plot_bin == TRUE) {
      x$k       <- x$center_k
      x$A       <- x$theta
      x$upper_A <- x$upper_bin
      x$lower_A <- x$lower_bin
      x$var_logA <- x$var_logbin
  } 
  dots <- function(...) {
      list(...)
  }

  
  additional_para <- dots(...)
  #print(additional_para)
  
  gg_color_hue = function( n ) {
    hues = seq( 15 , 375 , length = n + 1 )
    hcl( h = hues , l = 65 , c = 100 )[ 1 : n ]
  }
  cols = c( "grey25" , gg_color_hue( n = 3 ) )
  gray25 = cols[ 1 ]; red = cols[ 2 ]; green = cols[ 3 ]; blue = cols[ 4 ]
  
  if ("A" == plot[1]) {
      if (!is.null(high_deg_A))
          non_zero <- which(x$A > 10^-20 & x$k >= high_deg_A & x$k > 0)
      else 
          non_zero <- which(x$A > 10^-20 & x$k >= x$deg_threshold & x$k > 0) 
     
      if (!is.null(high_deg_A)) {
          v                   <- which(x$k[non_zero] == high_deg_A)  
          if (length(v) > 0) {
              new_var_log         <- x$var_logA[non_zero] #* (x$A[non_zero][v[1]]) ^ 2;
              x$A[non_zero]       <- x$A[non_zero] / x$A[non_zero][v[1]];    
              #print(x$var_logA[non_zero])
              #print(non_zero)
              #print(v)
              x$lower_A[non_zero] <- exp(log(x$A[non_zero]) - 2 * sqrt(new_var_log));
              x$upper_A[non_zero] <- exp(log(x$A[non_zero]) + 2 * sqrt(new_var_log));
              non_zero            <- x$A > 10^-20 & x$k >= high_deg_A &
                                   (!is.na(x$upper_A)) & (!is.infinite(x$upper_A))
          }
      } 
     if ((!is.null(max_A)) && (!is.null(min_A)))
          limit <- c(min(min_A,x$lower_A[non_zero]) , max(max_A,x$upper_A[non_zero]))
      else if (!is.null(min_A)) {
          limit <- c(min(min_A,x$lower_A[non_zero]) , max(x$upper_A[non_zero])) 
      }
      else if (!is.null(max_A)) {
          limit <- c(min(x$lower_A[non_zero]) , max(max_A,x$upper_A[non_zero]))
      }
      else { 
          limit <- c(min(x$lower_A[non_zero]) , max(x$upper_A[non_zero]))
      }
      
      if (is.null(additional_para$ylim))
        ylim <- limit  
      else ylim <- additional_para$ylim
   
      
      if (is.null(additional_para$xlim))
          xlim <- c(min(x$k[non_zero]) , max(x$k[non_zero]))  
      else xlim <- additional_para$xlim
      #print(xlim)
      
      temp   <- names(additional_para)
      ok_vec <- which(temp != "xlim" & temp != "ylim")
      #print(additional_para)
      additional_para_list <- ""
      for (i in ok_vec) {
          if (!is.character(additional_para[[i]]))  {
              additional_para_list <- paste0(additional_para_list ,temp[i]," = ",additional_para[[i]],",");
          } else  additional_para_list <- paste0(additional_para_list ,temp[i]," = '",additional_para[[i]],"',") 
      }
      final_para <- substr(additional_para_list,1,nchar(additional_para_list) - 1)
      #print(final_para)
      col_pa <- as.vector(col2rgb(col_point)) / 255
      
      eval(parse(text = paste('plot(x$k[non_zero][1], x$A[non_zero][1] , xlab = ifelse(!is.null(label_x),label_x, "Degree k"),
           ylab = ifelse(!is.null(label_y) , label_y,"Attachment function"),
           xlim = xlim , ylim = ylim , axes = FALSE, log = "xy" , 
           col = rgb(col_pa[1],col_pa[2],col_pa[3], shade_point), 
           mgp = c( 2.5 , 1 , 0 ),
           type = "n",', final_para, ')')))
      #magaxis(grid = TRUE, frame.plot = TRUE)
      eval(parse(text = paste('magaxis(grid = FALSE, frame.plot = TRUE,',final_para,')')));
      #xtick = seq(from = xlim[1], to = xlim[2],5)
      #axis(side = 1, at = xtick, labels = NULL, xlim = xlim, log = "x")
      if (TRUE == confidence) {
        order    <- order(x$k[non_zero])
        upper_f  <- x$upper_A[non_zero][order]
        lower_f  <- x$lower_A[non_zero][order]
        unique_x <- unique((x$k[non_zero])[order])
        upper_u  <- rep(0,length(unique_x))
        lower_u  <- upper_u    
        for (jjj in 1:length(unique_x)) {
          uu           <- which((x$k[non_zero])[order] == unique_x[jjj])
          #print(uu)
          upper_u[jjj] <- max(upper_f[uu])
          lower_u[jjj] <- min(lower_f[uu])
        }
        my_col <- as.vector(col2rgb(col_interval))
        my_col <- my_col / 255
        polygon(c(unique_x,rev(unique_x)) , c(upper_u,rev(lower_u)) , col = rgb(my_col[1],my_col[2],my_col[3],shade_interval),border = NA) 
        #arrows(x0 = x$k[non_zero] + 1, y0 = x$lower_A[non_zero], 
        #       x1 = x$k[non_zero] + 1, y1 = x$upper_A[non_zero], code = 3,angle = 90, 
        #       length = 0,
        #       col = rgb(0,0,0,shade_interval))
      }
      points(x$k[non_zero], x$A[non_zero] , pch = pch, col = rgb(col_pa[1],col_pa[2],col_pa[3], shade_point),...)
      if (TRUE == line) {
          alpha <- x$alpha
          #if (!is.null(names(x$loglinear_fit)))
          #    beta <-  x$loglinear_fit$coefficients[1]
          #else {
              index_one <- which(x$A[non_zero] == 1 & x$k[non_zero] != 0)[1] 
              beta      <- -alpha* log(x$k[non_zero][index_one])
              #print(beta)
          #}
          lines(x$k[non_zero], exp(beta) * (x$k[non_zero]) ^ alpha, lwd = 2, col = blue)
          
          #lines(c(1,x$k[non_zero] + 1),c(x$PA_offset , exp(beta) * (x$k[non_zero]) ^ alpha), lwd = 2, col = green)
      }
      

  }
  else if ("f" == plot[1]) {
      
      non_zero <- x$lower_f > 10^-20 & net_stat$increase > 0
      if (length(non_zero) <= 0)
        stop("There is no data. Please decrease high_deg") 
      # Plot the density of node fitnesses
      
      f_non <- x$f[non_zero]
      d     <- density(f_non)
      ok_d  <- d$x > 0
      max_x <- max(d$x[ok_d])
      min_x <- min(d$x[ok_d])
      
      #layout(cbind(1,2), width = c(4,1))
      plot(d$x[ok_d] , d$y[ok_d], col  = 2 , log = "x", lwd = 0, main = "", 
           xlab = "Fitness", ylab = "Density",xlim = c(min_x,max_x),
           mgp = c( 2.5 , 1 , 0 ), 
           axes = FALSE,
           pch = "",...)
      axis(side = 1, at=c(format(round(min_x, 2), nsmall = 2),
                          format(round(1, 2), nsmall = 2),
                          format(round(max_x, 2), nsmall = 2)),
           tcl = 0.5, ...)
      magaxis(side = 2,grid = FALSE, frame.plot = TRUE,
              logpretty = FALSE,...)
      
      #x <- c(format(min(f_non),digits = 1),1,5,10,15,format(max(f_non),digits = 4))
      u    <- smooth.spline(d$x, d$y, spar = 0.01)
      #polygon(d$x, d$y, col = red_fade, border=NA)
      ok_u <- u$x > 0 & u$y > 0 
      lines(u$x[ok_u], u$y[ok_u], col = "grey50",lwd = 2.5);
      #axis(3, at = 1,labels = "Mean = 1", las = 0,cex.axis = 2)
      gray <- rgb(0,0,0,1)
      #abline(v = median(f_non), lty = 5, lwd = 1.5,col = green)
      #abline(v = quantile(f_non,0.99), lty = 4, lwd = 1.5, col = blue)
      #mtext(at = median(f_non), side = 3,  
      #      format(round(median(f_non), 2), nsmall = 2), cex = 1)
      #mtext(at = quantile(f_non,0.99), side = 3,
      #      format(round(quantile(f_non,0.99), 2), nsmall = 2), 
      #      cex = 1)
      
      
      #plot.new()
      #par(mar=c(0, 0, 0, 0))
      #legend(legend = c("99th percentile" , "Median"), 
      #       text.col = c(gray),
      #       x = "left", col = c(green,blue),lwd = 1.5,
      #       lty = c(5,4),cex = 1, bty = "n")
  }
  else if ("true_f" == plot[1]) {
          #names(true_f) <- net_stat$node_id 
          #true_f        <-  true_f[names(x$f)] 
          true_f1   <- length(true_f[net_stat$node_id][net_stat$f_position]) * true_f[net_stat$node_id][net_stat$f_position]/
                       sum(true_f[net_stat$node_id][net_stat$f_position])
          if (FALSE == is.null(high_deg_f)) {
              non_zero <- x$lower_f[net_stat$f_position] > 10^-20 & true_f1 > 10^-20 & net_stat$increase[net_stat$f_position] > high_deg_f
          } else
              non_zero <- x$lower_f[net_stat$f_position] > 10^-20 & true_f1 > 10^-20 
          if (length(non_zero) <= 0)
             stop("There is no net_stat. Please decrease high_deg")  
          #print(non_zero)
          b        <- lm(true_f1[non_zero] ~ 0 + x$f[net_stat$f_position][non_zero])$coefficients[1]
          upper_f  <-  exp(log(b*x$f[net_stat$f_position][non_zero]) + 2 * sqrt(x$var_f[net_stat$f_position][non_zero] / x$f[net_stat$f_position][non_zero] ^ 2))
          lower_f  <-  exp(log(b*x$f[net_stat$f_position][non_zero]) - 2 * sqrt(x$var_f[net_stat$f_position][non_zero] / x$f[net_stat$f_position][non_zero] ^ 2))
          
          #print(lower_f)
          #print("---")
          #print(upper_f)
          
          if (is.null(additional_para$xlim)) {
              if ((!is.null(f_min)) && (!is.null(f_max))) {
                  if (TRUE == confidence)  
                      xlim <- c(min(c(f_min,lower_f,true_f1[non_zero])) , max(c(f_max,upper_f,true_f1)))
                  else xlim <- c(min(c(f_min,true_f1[non_zero])) , max(c(f_max,true_f1)))
              }
              else  if (!is.null(f_max)) {
                  if (TRUE == confidence)    
                      xlim <- c(min(c(lower_f,true_f1[non_zero])) , max(c(f_max,upper_f,true_f1)))
                  else xlim <- c(min(c(true_f1[non_zero])) , max(c(f_max,true_f1)))
              }    
              else if (!is.null(f_min)) {
                  if (TRUE == confidence)  
                      xlim <- c(min(c(f_min,lower_f,true_f1[non_zero])) , max(c(upper_f,true_f1)))  
                  else xlim <- c(min(c(f_min,true_f1[non_zero])) , max(c(true_f1)))  
              }
              else {   
                  if (TRUE == confidence)  
                      xlim <- c(min(c(lower_f,true_f1[non_zero])) , max(c(upper_f,true_f1)))  
                  else  xlim <- c(min(c(true_f1[non_zero])) , max(c(true_f1)))  
              }
          }
          else {
              xlim <- additional_para$xlim  
          }
          #print(xlim)  
          if (is.null(additional_para$ylim)) 
              ylim <- xlim
          else ylim <- additional_para$ylim
          temp   <- names(additional_para)
          ok_vec <- which(temp != "xlim" & temp != "ylim")
          #print(additional_para)
          additional_para_list <- ""
          for (i in ok_vec) {
            if (!is.character(additional_para[[i]]))  {
              additional_para_list <- paste0(additional_para_list ,temp[i]," = ",additional_para[[i]],",");
            } else  additional_para_list <- paste0(additional_para_list ,temp[i]," = '",additional_para[[i]],"',") 
          }
          final_para <- substr(additional_para_list,1,nchar(additional_para_list) - 1)
          
          #print(final_para)
          
          eval(parse(text = paste('plot(true_f1[non_zero][1], b * x$f[net_stat$f_position][non_zero][1] , xlim= xlim , 
                                   ylim = ylim,
                ylab = "Estimated fitness" , xlab = "True fitness" , log = "xy" , type = "n", axes = FALSE,
                                  ,
                 	mgp = c( 2.5 , 1 , 0 ) ,', 
                                  final_para, ')')))
          #magaxis(grid = TRUE, frame.plot = TRUE)
          eval(parse(text = paste('magaxis(grid = FALSE, frame.plot = TRUE,',final_para,')')));
          
         
          if (TRUE == confidence) {
              #x_point <- c(true_f1[non_zero],rev(true_f1[non_zero]))
              #y_point <- c(upper_f, rev(lower_f))
              order    <- order(true_f1[non_zero])
              upper_f  <- upper_f[order]
              lower_f  <- lower_f[order]
              unique_x <- unique(true_f1[non_zero][order])
              upper_u  <- rep(0,length(unique_x))
              lower_u  <- upper_u  
              for (jjj in 1:length(unique_x)) {
                  uu           <- which(true_f1[non_zero][order] == unique_x[jjj])
                  #print(uu)
                  upper_u[jjj] <- max(upper_f[uu])
                  lower_u[jjj] <- min(lower_f[uu])
              }
              my_col <- as.vector(col2rgb(col_interval))
              my_col <- my_col / 255
              polygon(c(unique_x , rev(unique_x)) , c(upper_u , rev(lower_u)) , col = rgb(my_col[1] , my_col[2] , my_col[3] , alpha = shade_interval) , border = NA)  
              #arrows(y0 = lower_f, x0 = true_f1[non_zero], y1 = upper_f, x1 = true_f1[non_zero], code = 3,angle = 90, length = 0,
              #       col = rgb(0,0,0,shade_interval))
          }
          abline(a=0,b = 1)
          if (FALSE == plot_true_degree) {
              col_fit <- as.vector(col2rgb(col_point)) / 255
              points(true_f1[non_zero],b * x$f[net_stat$f_position][non_zero] , pch = pch , col = rgb(col_fit[1],col_fit[2],
                                                                                                      col_fit[3],
                                                                                                      shade_point),...) 
          }
          else {  
              points(b * x$f[net_stat$f_position][non_zero],true_f1[non_zero] , pch = "" , ...)
              col         <- brewer.pal(9 , "Greens")
              order       <- order(net_stat$increase[net_stat$f_position][non_zero])
              col_seq     <- seq(min(net_stat$increase[net_stat$f_position][non_zero]) ,
                                 max(net_stat$increase[net_stat$f_position][non_zero]) , 9)
              a           <- sapply(net_stat$increase[net_stat$f_position][non_zero] , function(x)which(col_seq >= x)[1])
              a[is.na(a)] <- 9
              text(b * x$f[net_stat$f_position][non_zero] , true_f1[non_zero] , net_stat$increase[net_stat$f_position][non_zero] , col = col[a] , ...) 
          }
      }
}