# Packages needed

require("snowfall")
require("snowfall")
require("ggplot2")

##################################################################
###  With Rcpp package 

calcPar <- function(x, K, n.simu, r, type, path) {
library("Rcpp")  
sourceCpp(paste(path,"/codesCplusplus.cpp", sep = ""))
res <- f(K, n.simu, r, type)
return(res)
}

#######################################################
# parallel computing 
# initialization
choice.model <- c("IC", "IAC", "IACstar") 
K <- 51
B = 1000000000
nb.core <- 40
nb.simu.max.per.core <- 25000000
nb.jobs <- as.integer(B/(nb.simu.max.per.core))
B <- nb.jobs*nb.simu.max.per.core
path <- "/home/laurent/Documents/michel_lebreton/codes_R_simulations"

# the different sizes of population
taille.state = as.integer(10^c(0, 1, 2, 3, 4, 5, 6))

for (r in taille.state) { 

  sfInit(parallel = TRUE, cpus = nb.core)
  a = system.time(result.IC <- sfClusterApplyLB(1:nb.jobs, calcPar, K = K, n.simu = nb.simu.max.per.core, r = r, type = 0, path = path))
  print(a)
  a = system.time(result.IAC <- sfClusterApplyLB(1:nb.jobs, calcPar, K = K, n.simu = nb.simu.max.per.core, r = r, type = 1, path = path))
  print(a)
  a = system.time(result.IACstar <- sfClusterApplyLB(1:nb.jobs, calcPar, K = K, n.simu = nb.simu.max.per.core, r = r, type = 2, path = path))
  print(a)
  sfStop()

 # creation of the directory 
  if (!dir.exists(paste(path, "/r_", r, sep=""))) {
   dir.create(paste(path, "/r_", r, sep=""))
  }

  setwd(paste(path, "/r_", r, sep="")) 

 # save the results 
  save(result.IC, result.IAC, result.IACstar, file=paste("result_r", r, ".RData", sep=""))

  for (prob.model in choice.model){
    # select the results 
    result <- switch(prob.model,
                  IC = result.IC,
                  IAC = result.IAC,
                  IACstar = result.IACstar)

    pct.inversion.M1 = mean.Ut.in.M1 = mean.Ut.out.M1 = mean.Ut.M1 = 
      mean.pivot.in.M1 = mean.pivot.out.M1 = mean.pivot.M1 = 
      pct.inversion.M2 = mean.Ut.in.M2 = mean.Ut.out.M2 = mean.Ut.M2 = 
      mean.pivot.in.M2 = mean.pivot.out.M2 = mean.pivot.M2 = 
      mean.pivot.out.DO.M2 = mean.pivot.out.IO.M2 = mean.pivot.out.BD.M2 = 
      mean.pivot.out.BI.M2 = mean.pivot.out.B.M2 = numeric(K+1)
  
    for (k in 1:nb.jobs) {
      pct.inversion.M1 <- pct.inversion.M1 + result[[k]][,1]
      mean.Ut.in.M1 <- mean.Ut.in.M1 + result[[k]][,2] 
      mean.Ut.out.M1 <- mean.Ut.out.M1 + result[[k]][,3]
      mean.Ut.M1 <- mean.Ut.M1 + result[[k]][,4]  
      mean.pivot.in.M1 <- mean.pivot.in.M1 + result[[k]][,5]
      mean.pivot.out.M1 <- mean.pivot.out.M1 + result[[k]][,6]
      mean.pivot.M1 <- mean.pivot.M1 + result[[k]][,7]
  
      pct.inversion.M2 <- pct.inversion.M2 + result[[k]][,8]
      mean.Ut.in.M2 <- mean.Ut.in.M2 + result[[k]][,9]
      mean.Ut.out.M2 <- mean.Ut.out.M2 + result[[k]][,10]
      mean.Ut.M2 <- mean.Ut.M2 + result[[k]][,11]
      mean.pivot.in.M2 <-  mean.pivot.in.M2 + result[[k]][,12]
      mean.pivot.out.M2 <- mean.pivot.out.M2 + result[[k]][,13]
      mean.pivot.M2 <- mean.pivot.M2 + result[[k]][,14]
  
      mean.pivot.out.DO.M2 <- mean.pivot.out.DO.M2 + result[[k]][,15]
      mean.pivot.out.IO.M2 <- mean.pivot.out.IO.M2 + result[[k]][,16]
      mean.pivot.out.BD.M2 <- mean.pivot.out.BD.M2 + result[[k]][,17]
      mean.pivot.out.BI.M2 <- mean.pivot.out.BI.M2 + result[[k]][,18]
      mean.pivot.out.B.M2 <- mean.pivot.out.B.M2 + result[[k]][,19]
     }
  
    # we divide only to get the sharing in M2 case (pivot outside) 
    mean.pivot.out.DO.M2 <- mean.pivot.out.DO.M2/mean.pivot.out.M2
    mean.pivot.out.IO.M2 <- mean.pivot.out.IO.M2/mean.pivot.out.M2
    mean.pivot.out.BD.M2 <- mean.pivot.out.BD.M2/mean.pivot.out.M2
    mean.pivot.out.BI.M2 <- mean.pivot.out.BI.M2/mean.pivot.out.M2
    mean.pivot.out.B.M2 <- mean.pivot.out.B.M2/mean.pivot.out.M2
    
    # we divide the results by the number of simulations:
    pct.inversion.M1 <- pct.inversion.M1/B
    mean.Ut.in.M1 <- mean.Ut.in.M1/B
    mean.Ut.out.M1 <- mean.Ut.out.M1/B
    mean.Ut.M1 <- mean.Ut.M1/B
    mean.pivot.in.M1 <- mean.pivot.in.M1/B
    mean.pivot.out.M1 <- mean.pivot.out.M1/B
    mean.pivot.M1 <- mean.pivot.M1/B

    pct.inversion.M2 <- pct.inversion.M2/B 
    mean.Ut.in.M2 <- mean.Ut.in.M2/B
    mean.Ut.out.M2 <- mean.Ut.out.M2/B
    mean.Ut.M2 <- mean.Ut.M2/B
    mean.pivot.in.M2 <- mean.pivot.in.M2/B
    mean.pivot.out.M2 <- mean.pivot.out.M2/B
    mean.pivot.M2 <- mean.pivot.M2/B

    # we divide the results by the correct population
    taille_pop <- K * (2*r + 1) 
    taille_in <- (0:K)*(2*r + 1)
    taille_out <- taille_pop - taille_in
  
    mean.Ut.M1 <- mean.Ut.M1/taille_pop
    mean.pivot.M1 <- mean.pivot.M1/taille_pop
    mean.Ut.M2 <- mean.Ut.M2/taille_pop
    mean.pivot.M2 <- mean.pivot.M2/taille_pop 
  
    mean.pivot.in.M1 <- mean.pivot.in.M1/taille_in
    mean.pivot.in.M2 <- mean.pivot.in.M2/taille_in
    mean.Ut.in.M1 <- mean.Ut.in.M1/taille_in
    mean.Ut.in.M2 <- mean.Ut.in.M2/taille_in
  
    mean.pivot.out.M1 <- mean.pivot.out.M1/taille_out
    mean.pivot.out.M2 <- mean.pivot.out.M2/taille_out
    mean.Ut.out.M1 <- mean.Ut.out.M1/taille_out
    mean.Ut.out.M2 <- mean.Ut.out.M2/taille_out  
  
    # we correct the values for D=0 
    mean.pivot.M1[1] <- mean.pivot.out.M1[1]
    mean.pivot.M2[1] <- mean.pivot.out.M2[1]
    # we correct the values for D=K
    mean.pivot.M1[K+1] <- mean.pivot.in.M1[K+1]
    mean.pivot.M2[K+1] <- mean.pivot.in.M2[K+1]
  
    ##### creation of the files of results  
    if (!dir.exists(paste(getwd(), "/Figures_", prob.model, sep = ""))) {
      dir.create(paste(getwd(), "/Figures_", prob.model, sep = ""))
    }
  
    ###########
    # Inversions under M1 
    pdf(file = paste("Figures_", prob.model, "/inversion_", prob.model, "_M1.pdf", sep=""), 
      width = 10, height = 6)
  
    op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
    plot(0:K, pct.inversion.M1, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Probability of an election inversions under ", prob.model, " (Mechanism 1)", sep=""), 
       lty = 1, lwd = 3, col = "magenta", xaxt="n", las=2, 
       cex.axis = 0.8)
    axis(1)
    par(tcl =  -0.2)
    axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
    grid()
    legend("topright", legend = c("Maximum", "Back to initial probability"), lwd = 3, 
         pch = c(2, 6))
  
    # the max
    points(which.max(pct.inversion.M1)-1, max(pct.inversion.M1), col = "magenta", pch=2, cex=2)
    segments(which.max(pct.inversion.M1)-1, -0.5, which.max(pct.inversion.M1)-1, max(pct.inversion.M1), 
           lty = 2, col="grey")
    # the return to the beginning horizontal
    segments(which(pct.inversion.M1<pct.inversion.M1[1])[2] - 1, -0.5, 
           which(pct.inversion.M1<pct.inversion.M1[1])[2] - 1, 
           pct.inversion.M1[which(pct.inversion.M1<pct.inversion.M1[1])[2]], lty = 2, col="grey")
    points(which(pct.inversion.M1<pct.inversion.M1[1])[2] - 1, 
         pct.inversion.M1[which(pct.inversion.M1<pct.inversion.M1[1])[2]], 
         col = "magenta", pch=6, cex=2)
    par(op)
    dev.off()
  
    # Inversions under M2 
    pdf(file=paste("Figures_", prob.model, "/inversion_", prob.model, "_M2.pdf", sep=""), 
      width = 10, height = 6)
    op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
    plot(0:K, pct.inversion.M2, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Probability of an election inversions under ", prob.model, " (Mechanism 2)", sep=""), 
       lty = 1, lwd = 3, col = "magenta", xaxt="n", las=2, 
       cex.axis = 0.8)
    axis(1)
    par(tcl = -0.2)
    axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
    grid()
    dev.off()
  
    #############
    # utilities inside/outside the coalition under M1 
    pdf(file=paste("Figures_", prob.model, "/utilities_", prob.model, "_M1.pdf", sep=""), 
      width = 10, height = 6)
    op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
    plot(0:K, mean.Ut.M1, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Utility under ", prob.model, " (Mechanism 1)", sep=""), 
       lty = 1, lwd = 3, col = "cyan", xaxt = "n", las = 2, 
       ylim = range(c(mean.Ut.in.M1, mean.Ut.out.M1), na.rm = T), 
       cex.axis = 0.8)
    axis(1)
    par(tcl =  -0.2)
    axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
    grid()
    lines(0:K, mean.Ut.out.M1, type = "l", lty = 3, lwd = 3, col = "cyan")
    lines(0:K, mean.Ut.in.M1, type = "l", lty = 2, lwd = 3, col = "cyan")
  
    legend("topright", legend = c("Minimum", "Back to initial level", "Total", "Inside", "Outside"), 
         lwd = 3, lty = c(NA, NA, 1:3),
         pch = c(6, 2, NA, NA, NA))
  
    # the min of the total 
    points(which.min(mean.Ut.M1)-1, min(mean.Ut.M1, na.rm=T), col = "cyan", pch=6, cex=2)
    segments(which.min(mean.Ut.M1)-1, -0.5, which.min(mean.Ut.M1)-1, min(mean.Ut.M1, na.rm=T), 
           lty = 2, col="grey")
    # the return to the beginning horizontal
    segments(which(mean.Ut.M1>mean.Ut.M1[1])[1] - 1, -0.5, which(mean.Ut.M1>mean.Ut.M1[1])[1] - 1, 
           mean.Ut.M1[which(mean.Ut.M1>mean.Ut.M1[1])[1]], lty = 2, col="grey")
    points(which(mean.Ut.M1>mean.Ut.M1[1])[1] - 1, 
         mean.Ut.M1[which(mean.Ut.M1>mean.Ut.M1[1])[1]], col = "cyan", pch=2, cex=2)
    dev.off()
  
    #############
    # utilities inside/outside the coalition under M2 
    pdf(file=paste("Figures_", prob.model, "/utilities_", prob.model, "_M2.pdf", sep=""), 
      width = 10, height = 6)
    op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
    plot(0:K, mean.Ut.M2, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Utility under ", prob.model, " (Mechanism 2)", sep=""), 
       lty = 1, lwd = 3, col = "cyan", xaxt = "n", las = 2, 
       ylim = range(c(mean.Ut.in.M2, mean.Ut.out.M2), na.rm = T), cex.axis = 0.8)
    axis(1)
    par(tcl =  -0.2)
    axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
    grid()
    lines(0:K, mean.Ut.out.M2, type = "l", lty = 3, lwd = 3, col = "cyan")
    lines(0:K, mean.Ut.in.M2, type = "l", lty = 2, lwd = 3, col = "cyan")
  
    legend("topright", legend = c("Total", "Inside", "Outside"), lwd = 3, lty = 1:3)
  
    # the min of the total 
    dev.off()
  
    #############
    # pivot inside/outside the coalition under M1
    pdf(file=paste("Figures_", prob.model, "/pivot_", prob.model, "_M1.pdf", sep=""), 
      width = 12, height = 8)
    op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
    plot(0:K, mean.pivot.in.M1, type = "l", lty = 2, lwd = 3, col = "green", 
       ylim = range(c(mean.pivot.M1, mean.pivot.in.M1, mean.pivot.out.M1), na.rm = T),
       ylab = paste("Probability to be pivot ", " (Mechanism 1)", sep=""), 
       xlab = "Number of states inside the coalition", las=2, 
       xaxt = "n", cex.axis = 0.8)
    axis(1)
    par(tcl =  -0.2)
    axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
    grid()
    lines(0:K, mean.pivot.out.M1, type = "l", lty = 3, lwd = 3, col = "green")
    lines(0:K, mean.pivot.M1, type = "l", lty = 1, lwd = 3, col = "green")
    legend("topright", legend = c("Total", "Inside", "Outside"), lwd = 3, lty = 1:3)
    par(op)
    dev.off()
  
    #############
    # pivot inside/outside the coalition under M2
    pdf(file=paste("Figures_", prob.model, "/pivot_", prob.model, "_M2.pdf", sep=""), 
      width = 12, height = 8)
    op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
    plot(0:K, mean.pivot.in.M2, type = "l", lty = 2, lwd = 3, col = "green", 
       ylim = range(c(mean.pivot.M2, mean.pivot.in.M2, mean.pivot.out.M2), na.rm = T),
       ylab = paste("Probability to be pivot ", " (Mechanism 2)", sep=""), 
       xlab = "Number of states inside the coalition", las=2, 
       xaxt = "n", cex.axis = 0.8)
    axis(1)
    par(tcl =  -0.2)
    axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
    grid()
    lines(0:K, mean.pivot.out.M2, type = "l", lty = 3, lwd = 3, col = "green")
    lines(0:K, mean.pivot.M2, type = "l", lty = 1, lwd = 3, col = "green")
    legend("topright", legend = c("Total", "Inside", "Outside"), lwd = 3, lty = 1:3)
    par(op)
    dev.off()
  
    #############
    # pivot outside the coalition under M2
    data.ggplot <- data.frame(Sharing = c(mean.pivot.out.DO.M2[-(K+1)], 
                                        mean.pivot.out.IO.M2[-(K+1)], 
                                        mean.pivot.out.BD.M2[-(K+1)], 
                                        mean.pivot.out.BI.M2[-(K+1)], 
                                        mean.pivot.out.B.M2[-(K+1)]),
                            Sector = c(rep("DO", K), rep("IO", K), rep("BD", K), 
                                       rep("BI", K), rep("B", K)),
                            L = rep(0:(K-1), 5))
  
    data.ggplot$Sharing[is.na(data.ggplot$Sharing)] <- 0
    myplot <- ggplot(data.ggplot, aes(x = L, y = Sharing, fill = Sector)) + 
      geom_area(alpha = 0.6 , size = 1, colour = "black")
  
    ggsave(filename=paste("Figures_", prob.model, "/pivot_", prob.model, "_outside_M2.pdf", sep=""), 
         plot = myplot)
   }
}