#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#Here I implement a Grid search method for finding all LQRE in a two by two strategic form game:
#
#Mark Fey
#4/9/2009
#
#Adapted from code by
#Jeremy Kedziora
#3/5/2009
#In these boxes the numbers list the order that the utilities appear in the vector ui
#   L     R
#  _____________
#U |  1  |  2  |
#  |     |     |
#  -------------
#D |  3  |  4  |
#  |     |     |
#  -------------

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
rm(list=ls())




#@@@@@@@@@@@@@@@@@@@@@@@
#input the utilties here
#@@@@@@@@@@@@@@@@@@@@@@@
u1<-c(2,0,0,1)
u2<-c(2,0,0,1)

# Battle of Sexes
#u1<-c(3,0,0,1)
#u2<-c(1,0,0,2)

#u1<-c(2,2,0,3)
#u2<-c(2,2,0,1)


#@@@@@@@@@@@@@@@@@@@@@@@
#other values defined here.
#lambda values in plot from 0 to lambdamax, by step lambdainc
#grid size for equilibrium search is gridsize
#@@@@@@@@@@@@@@@@@@@@@@@
lambdamax<-15
lambdainc<-0.01
gridsize<-1000

#@@@@@@@@@@@@@@@@@@@@@@@@@
#Define the prob of choosing between alt A and alt B, given some parameter
#@@@@@@@@@@@@@@@@@@@@@@@@@
ProbA<-function(A,B,param){
    exp(param*A)/(exp(param*A)+exp(param*B))
}


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#Then the best response function for player 2 given player one's probability - we'll actually solve for p and then compute q directly
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
BR2<-function(p,u1,u2,lambda){
    ProbA((p*u2[1]+(1-p)*u2[3]),(p*u2[2]+(1-p)*u2[4]), lambda)
}



#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#First the best response function for player 1 given player two's probability
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
BR1<-function(p,u1,u2,lambda){
    q<-BR2(p,u1,u2,lambda)
    ProbA((q*u1[1]+(1-q)*u1[2]),(q*u1[3]+(1-q)*u1[4]), lambda)
}




#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#We'll compute the fixed point in p with this
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
R1<-function(p,u1,u2,lambda){
    BR1(p,u1,u2,lambda)-p   
}


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#Here for fixed lambda we'll compute the LQRE
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
LQRE<-function(u1,u2,lambda,npvals){
    
    #Here we create the grid for the p:
    pvals<-seq(0,1,length=npvals)
    
    #Here we bracket the roots
    dif<-(-1)*R1(pvals,u1,u2,lambda)
    dif.sign<-dif>=0
    change<-(dif.sign*1)[1:(npvals-1)]-(dif.sign*1)[2:(npvals)]
    
    #for any of these values, we need to take where a non-zero appears and then the next one:
    Lo<-subset(c(abs(change),0)*seq(1,npvals),c(abs(change),0)*seq(1,npvals)!=0)
    Hi<-Lo+1

    roots<-matrix(NA,nrow=length(Lo),ncol=2)
    for(i in 1:length(Lo)){
        roots[i,1]<-uniroot(R1,c(pvals[Lo[i]],pvals[Hi[i]]),u1=u1,u2=u2,lambda=lambda,tol=1e-15)$root
        roots[i,2]<-BR2(roots[i,1],u1,u2,lambda)
    }
    
    roots
    
}


#@@@@@@@@@@@@@@@@@
#Main Program
#@@@@@@@@@@@@@@@@@

#Generate the lambda values, starting at 0 and ending at lambdamax, step lambdainc
lambda<-seq(0,lambdamax,by=lambdainc)

#Initialize the matrix to hold the values that we will plot
QRE.ALL<-matrix(NA,nrow=1,ncol=3)

#Build the QRE.ALL matrix 
#The QRE function will return a matrix with one row for each QRE at the lambda value
#The first column is the lambda value and the second and third column are the p and q
# values of the equilibrium.
for(j in 1:length(lambda)){
    QRE.ALL<-rbind(QRE.ALL,cbind(lambda[j],LQRE(u1,u2,lambda[j],gridsize)))
}


#Plot the result.
#Create two graphs.  One is the p-value and one is the q-value.
split.screen(c(2,1))
screen(1)
plot("",ylim=c(0,1),xlim=c(0,lambdamax))
for(j in 1:length(QRE.ALL[,1])){
    points(lambda[j],QRE.ALL[j,2],pch=20,cex=0.1)
}
    

screen(2)
plot("",ylim=c(0,1),xlim=c(0,lambdamax))
for(j in 1:length(QRE.ALL[,1])){
    points(lambda[j],QRE.ALL[j,3],pch=20,col="grey",cex=0.1)
}