#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 23
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

# From my c23.txt:

# c23.txt
# 18 Apr 2013

# Kelly criterion, cont'd.


GamblersRuinS := proc(p, start, goal, strat) local i, amt, bet:
    if start < 0 or start > goal then:
        return FAIL:
    fi:
    amt := start:
    for i from 0 do:
        if amt = goal then:
            return true, i:
        fi:
        if amt = 0 then:
            return false, i:
        fi:
        bet := min(strat(amt), goal-amt):
        if bet > amt or bet <= 0 then:
            return FAIL, amt:
        fi:

        if LC(p) then:
            amt := amt + bet:
        else:
            amt := amt - bet:
        fi:
    od:
end:

GRSMonteCarlo := proc(p, start, goal, strat, N) local L, i, wins, losses, dur:
    NumericEventHandler(division_by_zero=proc(operator, operands, defVal) return infinity: end):
    L := [seq([GamblersRuinS(p, start, goal, strat)], i=1..N)]:
    wins :=  select(x->x[1], L):
    losses := select(x->not x[1], L):
    dur := add(i[2], i in L):
    return nops(wins)/N, add(i[2], i in wins)/nops(wins), add(i[2], i in losses)/nops(losses), dur/N:
end:

VWS := proc(p, goal, strat) local v, eqs, vars, solns:
    vars := [seq(v[i], i=0..goal)]:
    eqs := {v[0]=0, v[goal]=1}:
    eqs := eqs union {seq(p*v[i+min(strat(i), goal-i)] + (1-p)*v[i-min(strat(i), goal-i)] = v[i], i=1..goal-1)}:
    solns := solve(eqs, vars):
    return subs(convert(solns[1], set), vars):
end:

VDS := proc(p, goal, strat) local v, eqs, vars, solns:
    vars := [seq(v[i], i=0..goal)]:
    eqs := {v[0]=0, v[goal]=0}:
    eqs := eqs union {seq(p*v[i+min(strat(i), goal-i)] + (1-p)*v[i-min(strat(i), goal-i)] +1 = v[i], i=1..goal-1)}:
    solns := solve(eqs, vars):
    return subs(convert(solns[1], set), vars):
end:

Dominates := proc(L1, L2):
    return not evalb(true in {seq(evalb(L1[i] < L2[i]), i=1..nops(L1))}):
end:

AllIntegerFuncsLess := proc(f, n) local i:
    return map(x->(y->x[y]), AllListsLess([seq(f(i), i=1..n)])):
end:

AllListsLess := proc(L) local n,i,j: option remember:
    n := nops(L):
    if n = 0 then:
        return {[]}:
    else:
        return {seq(seq([op(i), j], i=AllListsLess([op(1..n-1, L)])), j=1..L[n])}:
    fi:
end:

AllStrats := proc(n):
    return AllIntegerFuncsLess(x->min(x, n-x), n-1):
end:

BestStrat := proc(p, n) local S, maximal, current, new, i:
    S := AllStrats(n):
    maximal := {S[1]}:
    current := VWS(p, n, S[1]):
    for i from 2 to nops(S) do:
        new := VWS(p,n,S[i]):
        if new = current then:
            maximal := maximal union {S[i]}:
        elif Dominates(new, current) then:
            maximal := {S[i]}:
            current := new:
        fi:
    od:
    return maximal, current:
end:

#############
# Problem 1 #
#############

TimidVsBold := proc(n, p) local i:
    return [seq({solve({VWS(p, n, x->1)[i+1] >= VWS(p, n, x->x)[i+1], 0 <= p, p <= 1}, {p})}, i=1..n-1)]:
end:

# solve choked when n >= 15, but for n <= 14 timid is better than bold
# whenever p >= 1/2. (Also, it doesn't matter at 0 and 1.)

#############
# Problem 2 #
#############

Neighbors := proc(n,f) local i,j, k:
    return map(x->(y->x[y]), select(IsLegal, {seq(seq([seq(f(j), j=1..i-1), f(i)+k, seq(f(j), j=i+1..n-1)], i=1..n-1), k={+1, -1})})):
end:

IsLegal := proc(L) local i:
    return not evalb(true in {seq(evalb(L[i] <= 0), i=1..nops(L)), seq(evalb(L[i] > min(i, nops(L)+1-i)), i=1..nops(L))}):
end:

BestNeighbor := proc(n,p,f,i) local best, val, new, newval:
    best := f:
    val := VWS(p, n, f)[i]:
    for new in Neighbors(n, f) do:
        newval := VWS(p, n, new)[i]:
        if newval > val then:
            best := new:
            val := newval:
        fi:
    od:
    return best:
end:

#############
# Problem 3 #
#############

BestStratFast := proc(p, N, i) local f, g, j:
    f := x->1:
    while true do:
        g := BestNeighbor(N, p, f, i):
        if g = f then:
            return [seq(f(j), j=1..N-1)]:
        else:
            f := g:
        fi:
    od:
    return [seq(f(j), j=1..N-1)]:
end:

#############
# Problem 4 #
#############

# It seems to work.