#Nathan Fox
#Homework 24
#I give permission for this work to be posted online

#Read procedures from class/last homework
read(`C24.txt`):

#Help procedure
Help:=proc():
    print(` GainLS(m, p, S) , SimGain(m, p, S, M) , Neigbors(S) `):
    print(` BestNeighbor(m, p, S, M) , BestStrat(m, p, M) `):
end:

##PROBLEM 1##

#GainLS(m, p, S): an abbeviated (faster) version of
#PlayLSterse(m, p, S) that only computes the final gain
#(not the duration)
GainLS:=proc(m, p, S) local m1, p1, ga, db:
    m1:=m:
    p1:=p:
    ga:=0:
    
    while m1+p1>0 do
        if p1 = 0 or m1-p1>=S[p1] then
            return ga:
        else
            db:=rand(1..p1+m1)():
            if db <= m1 then
                m1:=m1-1:
                ga:=ga-1:
            else
                p1:=p1-1:
                ga:=ga+1:
            fi:
        fi:
    od:
    return ga:
end:

#SimGain(m, p, S, M): runs a Shepp Urn game with m minus balls and
#p plus balls, M times, and computes the average.
SimGain:=proc(m, p, S, M) local i:
    return add(GainLS(m, p, S), i=1..M)/M:
end:

##PROBLEM 2##

#Neigbors(S): inputs a strategy, and outputs the set of strategies
#where only one component is changed, by either 1 or -1 (and the
#entries are never negative), so there are at most 2*nops(S)
#neighbors of a strategy S.
Neighbors:=proc(S) local ret, i:
    ret:={}:
    for i from 1 to nops(S) do
        if S[i] > 0 then
            ret:=ret union {[op(1..i-1, S), S[i]-1, op(i+1..nops(S), S)]}:
        fi:
        ret:=ret union {[op(1..i-1, S), S[i]+1, op(i+1..nops(S), S)]}:
    od:
    return ret:
end:

##PROBLEM 3##

#BestNeighbor(m, p, S, M): inputs positive integers m and p, a
#strategy, S, and a large integer M (about a 1000), and outputs,
#the neighboring strategy that gives the best score of
#SimGain(m, p, S', M) (among all the neighbors S' of S) if it is
#better than that of S, or returns S, if all the neighbors have a
#lower score.
BestNeighbor:=proc(m, p, S, M) local NS, champ, rec, N, score:
    NS:=Neighbors(S):
    champ:=S:
    rec:=SimGain(m, p, S, M):
    for N in NS do
        score:=SimGain(m, p, N, M):
        if score > rec then
            rec:=score:
            champ:=N:
        fi:
    od:
    return champ:
end:

##PROBLEM 4##

#BestStrat(m, p, M): starts with the "chicken" strategy S=[0$p]
#and by the above "hill-climbing" finds a locally best strategy
#for m minus balls and p plus balls.
BestStrat:=proc(m, p, M) local S, T:
    S:=[0$p]:
    while true do
        T:=BestNeighbor(m, p, S, M):
        if T = S then
            return S:
        fi:
        S:=T:
    od:
end:

#BestStrat(10,10,1000) ran for a REALLY long time, but eventually returned [1, 2, 3, 3, 4, 3, 3, 3, 2, 1]
#BestStrat(14,16,1000) ran for a REALLY REALLY long time, but eventually returned [1, 2, 2, 2, 3, 2, 4, 3, 3, 3, 3, 3, 3, 1, 1, 0]