rm(list=ls())
setwd("C:/Users/Yossi/Documents/Various R programs")
party=c("likud", "bluewhite", "meshutefet",
        "shas", "aguda", "yamina", "labor", 
        "liberman", "democrats", "otzma", 
        "noam", "others")
votes=c(1140370, 1125881, 337108, 
        258275, 249209, 283910, 190870, 
        173004, 135529, 133211, 
        75223, 33333)
k22hypo=data.frame(party, votes)
k22hypo

rm(party, votes)

combine_votes=function(current, party1, party2){ 
  # party1="likud"; party2="yamin_united"; current=suragg
  if((party1 %in% current$party) & (party2 %in% current$party)){
    nr=nrow(current)+1
    current[nr, 1]=paste(party1, party2)#, sep="+")
    w=which(current$party %in% c(party1, party2))
    current[nr, 2]=sum(current$votes[w])
    current[nr, 3]=sum(current$mandats[w])
    for (j in 2:1){
      current=current[-w[j],]
    }
  }
  return(current)
}


# surplus agreements
# liberman + bluewhite
# likud + yamina
# labor + democrats
# shas + aguda
# noam + otzma


surplus_aggregate=function(knstsurp){ # knstsurp=knst
  suragg=knstsurp
  suragg$party=as.character(suragg$party)
  suragg=combine_votes(suragg, "liberman", "bluewhite")
  suragg=combine_votes(suragg, "likud", "yamina")
  suragg=combine_votes(suragg, "labor", "democrats")
  suragg=combine_votes(suragg, "shas", "aguda")
  suragg=combine_votes(suragg, "noam", "otzma")
  rownames(suragg)=1:nrow(suragg)
  return(suragg)
}


add_mandat=function(add, w){ # add=knst
  add=add[w,]
  add$pay=add$votes/(add$mandats+1)
  w1=which(add$pay==max(add$pay))
  if(length(w1)>1){
    w1=rbinom(1,1,0.5)+1
  }
  select=as.character(add[w1, "party"])
  return(select)
}

bader_ofer=function(knst){ # knst=knesset
  bo=surplus_aggregate(knst) # bo=suragg
  bo$pay=bo$votes/(1+bo$mandats)
  w=which(bo$pay==max(bo$pay))
  parties=strsplit(bo$party[w], split=" ")[[1]]
  w=which(knst$party %in% parties)
  select=add_mandat(knst,w)
  w=which(knst$party==select)
  knst$mandats[w]=knst$mandats[w]+1
  return(knst)
}

the_knesset=function(results){ # results=thepoll
  pct=sum(results$votes)*0.0325
  w=which(results$votes>pct)
  results=results[w, ]
  knesset=results
  
  # calculate moded
  moded=sum(results$votes)/120
  
  # assign mandats by moded
  knesset$mandats=floor(results$votes/moded)
  remained=120-sum(knesset$mandats)
  if(remained<0){
    print("negative number of mandats remained")
    break
  }
  
  for (j in 1:remained){
    knesset=bader_ofer(knesset)
  }
  return(knesset)  
}

population=k22hypo
population=population[order(population$party),]
population$start=1
population$end=cumsum(population$votes)
population$start[2:nrow(population)]=population$end[1:(nrow(population)-1)]+1
population$party=as.character(population$party)
pop_size=sum(population$votes)
n_parties=nrow(population)
population


real_res=the_knesset(k22hypo)
w=which(!(k22hypo$party %in% real_res$party))
did_not_pass=k22hypo[w,]
did_not_pass$mandats=0
real_res=rbind(real_res, did_not_pass)
real_res=real_res[-n_parties,]
real_res=real_res[order(real_res$party),]

rm(k22hypo, did_not_pass, w)


make_sim_res=function(){
  res=as.data.frame(matrix(0, nrow=1, ncol=6))
  names(res)=c("s", "mean_diff",  "max_diff", 
               "med_diff", "pass_error", "did_not_pass_error")
  return(res)
}


cat(rep("\n", 50))
sample_size=c(500, 1000, 2000, 4000, 8000)
nsims=2000

sim_final_res=make_sim_res()
names(sim_final_res)[1]="sample_size"

set.seed(1310)
for(s in 1:length(sample_size)){
  sim_res=make_sim_res()
  for (j in 1:nsims){
    if(j %% 100 == 0){
      string=paste("sample size=", sample_size[s], "  j=", j, sep="")
      print(string)
    }
    the_sample=sample(1:pop_size, sample_size[s])  
    the_poll=" "
    for (k in 1:length(the_sample)){
      selected=rep(the_sample[k], n_parties)
      w=which(selected>=population$start & selected<=population$end)
      the_poll[k]=population$party[w]
    }
    the_poll=data.frame(table(the_poll))
    names(the_poll)=c("party", "votes")
    poll_mandats=the_knesset(the_poll)[,-2]
    names(poll_mandats)[2]="poll_res"
    w=which(poll_mandats$party=="other")
    poll_mandats=merge(real_res, poll_mandats, by="party", all.x=TRUE)
    w=which(is.na(poll_mandats$poll_res))
    if(length(w)>0){
      poll_mandats$poll_res[w]=0
    }
    diff=abs(poll_mandats$poll_res-poll_mandats$mandats)
    sim_res[j,1]=j
    sim_res$mean_diff[j]=mean(diff)
    sim_res$max_diff[j]=max(diff)
    sim_res$med_diff[j]=median(diff)
    sim_res$pass_error[j]=length(which(poll_mandats$mandats==0 & poll_mandats$poll_res>0))
    sim_res$did_not_pass_error[j]=length(which(poll_mandats$mandats>0 & poll_mandats$poll_res==0))
  }
  sim_final_res[s,1]=sample_size[s]
  sim_final_res[s,2:6]=apply(sim_res[, 2:6], 2, mean)
}