entry-game.R

library(data.table)
library(ggplot2)

rm(list=ls())
#set.seed(1)

# solver parameters
eps <- 1e-2     # need an epsilon for the root finder
mingap <- 0.3

# hyperparameters
a <- 1
b <- 1
sigma <- 1

# functions
u <- function(m) a*m - exp(-b*m)
ct <- function(delta) -1/delta
dct <- function(delta) 1/(delta^2)

Edu <- function(delta,W0) a + b*exp(-b*(W0 + ct(delta)) + 0.5*b^2*sigma^2*delta)
Ed2u <- function(delta,W0) -b^2 * exp(-b*(W0 + ct(delta)) + 0.5*b^2*sigma^2*delta)

Euopen <- function(delta,W0) a*(W0+ct(delta))-exp(-b*(W0+c(delta))+0.5*b^2*sigma^2*delta)
Eubrid <- function(delta,xl,xh,Wl,Wh) {
    dbar <- xh - xl - delta
    arg <- delta*(Wh-Wl)/(xh-xl) + Wl + ct(delta) + ct(dbar)
    a*arg - exp(-b*arg + 0.5*b^2*sigma^2*(delta*dbar)/(xh-xl))
}

opencrit <- function(delta,W0) (2*dct(delta))/(sigma^2) + Ed2u(delta,W0)/Edu(delta,W0)
bridcrit <- function(delta,xl,xh,Wl,Wh) {
    dbar <- xh - xl - delta
    (2*(Wh-Wl)/(sigma^2 * (dbar-delta))) + 2*(dbar-delta)/(sigma^2*(xh-xl))* (dct(delta)-dct(dbar)) + Ed2u(delta,Wh)/Edu(delta,Wh)
}

openjump <- function(delta) sigma*sqrt(delta)*rnorm(1)
bridjump <- function(delta,xl,xr,yl,yr) yl + delta*(yr-yl)/(xr-xl) + sigma*sqrt((delta*(xr-xl-delta))/(xr-xl))*rnorm(1)

dt <- data.table(m=seq(-5e-1,5e-1,length.out=1e4))
dt[, delta := seq(5e-1,1e2,length.out=.N)]
dt[, u := u(m)]
dt[, Edu := Edu(delta,0)]
dt[, Ed2u := Ed2u(delta,0)]
dt[, R := Ed2u/Edu]
dt[, oc := opencrit(delta,20)]

dt[, db := seq(0.5+eps,1-eps,length.out=.N)]
dt[, bc := bridcrit(db,0,1,0,1)]

ggplot(dt,aes(db,bc)) + geom_point()

#break()


xs <- c(0)
Ws <- c(0)
jumps <- NULL
gaps <- NULL

# setup
n <- 1e3
info <- data.table(id=1:(n+2),xl=-Inf,xr=Inf,
                   Wl=as.numeric(NA),Wr=as.numeric(NA))
d0 <- uniroot(opencrit,c(5e-2,1e2),W0=0)$root
Eu0 <- Euopen(d0,0)
info[1, `:=`(xr=0,Wr=0,delta=d0,Eu=Eu0,xd=-d0)]
info[2, `:=`(xl=0,Wl=0,delta=d0,Eu=Eu0,xd=d0)]

for (i in 3:(n+2)) {
    #print(i)
    #print(info)
    plot(c(info[is.finite(xr),xr],info[!is.na(xd),xd]),
         c(info[is.finite(xr),Wr],info[!is.na(xd),Eu]),
         pch=c(rep(15,i-2),rep(20,i-1)),ylim=c(-15,50))

    # pick best option, and split
    bestidx <- info[,order(Eu,decreasing=TRUE)[1]]
    if(is.infinite(info[bestidx,xl])) {
        Wnew <- info[bestidx,Wr] + openjump(info[bestidx,delta])
    } else if (is.infinite(info[bestidx,xr])) {
        Wnew <- info[bestidx,Wl] + openjump(info[bestidx,delta])
    } else {
        Wnew <- bridjump(info[bestidx,delta],
                         info[bestidx,xl],info[bestidx,xr],
                         info[bestidx,Wl],info[bestidx,Wr])
    }

    # split old interval into left part (bestidx) and right part (i)
    info[i, `:=`(xl=info[bestidx,xd],xr=info[bestidx,xr],
                 Wl=Wnew,Wr=info[bestidx,Wr])]
    info[bestidx, `:=`(xr=xd,Wr=Wnew)]

    # get new interval expectations for the lefter interval (bestidx)
    if (is.infinite(info[bestidx,xl])) {
        d0 <- uniroot(opencrit,c(5e-2,1e2),W0=info[bestidx,Wr])$root
        Eu0 <- Euopen(d0,info[bestidx,Wr])
        info[bestidx, `:=`(delta=d0,Eu=Eu0,xd=xr-d0)]
    } else {
        if(info[bestidx,xr-xl>mingap]){
            if(info[bestidx,Wr>Wl]) {
                rootl <- info[bestidx,(xr-xl)/2 + eps]
                rootr <- info[bestidx,(xr-xl) - eps]
            } else {
                rootl <- info[bestidx,eps]
                rootr <- info[bestidx,(xr-xl)/2 - eps]
            }
            #print("beep")
            #print(info[bestidx,xr-xl])
            gaps <- c(gaps,(info[bestidx,xr-xl]))
            d0 <- uniroot(bridcrit,c(rootl,rootr),
                          xl=info[bestidx,xl],xh=info[bestidx,xr],
                          Wl=info[bestidx,Wl],Wh=info[bestidx,Wr])$root
            Eu0 <- Eubrid(d0,info[bestidx,xl],info[bestidx,xr],
                          info[bestidx,Wl],info[bestidx,Wr])
        } else {
            d0 <- info[bestidx,(xr-xl)/2]
            Eu0 <- -Inf
        }
        info[bestidx, `:=`(delta=d0,Eu=Eu0,xd=xl+d0)]
    }

    # get new interval expectations for the righter interval (i)
    if (is.infinite(info[i,xr])) {
        d0 <- uniroot(opencrit,c(5e-2,1e2),W0=info[i,Wl])$root
        Eu0 <- Euopen(d0,info[i,Wl])
        info[i, `:=`(delta=d0,Eu=Eu0,xd=xl+d0)]
    } else {
        if(info[i,xr-xl>mingap]){
            if(info[i,Wr>Wl]) {
                rootl <- info[i,(xr-xl)/2 + eps]
                rootr <- info[i,(xr-xl) - eps]
            } else {
                rootl <- info[i,eps]
                rootr <- info[i,(xr-xl)/2 - eps]
            }
            #print("boop")
            #print(info[i,xr-xl])
            gaps <- c(gaps,(info[i,xr-xl]))
            d0 <- uniroot(bridcrit,c(rootl,rootr),
                          xl=info[i,xl],xh=info[i,xr],
                          Wl=info[i,Wl],Wh=info[i,Wr])$root
            Eu0 <- Eubrid(d0,info[i,xl],info[i,xr],
                          info[i,Wl],info[i,Wr])
        } else {
            d0 <- info[i,(xr-xl)/2]
            Eu0 <- -Inf
        }
        info[i, `:=`(delta=d0,Eu=Eu0,xd=xl+d0)]
    }
    Sys.sleep(0.01)
    #readline()
}

print(info)
#plot(c(info[is.finite(xr),xr],info[!is.na(xd),xd]),
     #c(info[is.finite(xr),Wr],info[!is.na(xd),Eu]),
     #pch=c(rep(19,n-1),rep(1,n)))


#ggplot(dt,aes(m,u)) + geom_point()
#ggplot(dt,aes(delta,Edu)) + geom_point()
#ggplot(dt,aes(delta,Ed2u)) + geom_point()
#ggplot(dt,aes(delta,R)) + geom_point()
#ggplot(dt,aes(delta,oc)) + geom_point()