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

# Nathaniel Shar
# HW 24
# Experimental Mathematics

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

# From my c24.txt:

with(RandomTools[MersenneTwister]):
randomize():
f2L := proc(f, n) local i:
    return [seq(f(i), i=1..n-1)]:
end:

L2f := proc(L):
    return x->L[x]:
end:

GetMax := proc(S, f:=(x->x)) local champion, record, new, newval, i:
    champion := S[1]:
    record := f(champion):
    for i from 2 to nops(S) do:
        new := S[i]:
        newval := f(new):
        if newval > record then:
            champion := new:
            record := newval:
        fi:
    od:
    return [champion, record]:
end:

# Probability of going from i dollars to N dollars by time T, with
# best strategy
BestSDD := proc(N, i, p, T) local j: option remember:
    if i = N then:
        return(true, 1):
    elif i = 0 then:
        return(false, 0):
    elif T = 0 then:
        return(false, 0):
    else:
        return GetMax({seq(j, j=1..min(i, N-i))}, x->(p*BestSDD(N, i+x, p, T-1)[2] + (1-p)*BestSDD(N, i-x, p, T-1)[2])):
    fi:
end:

SimulateOptimalGambling := proc(N, i, p, T) local cap, j, bet:
    cap := i:
    for j from T to 1 by -1 do:
        bet := BestSDD(N, cap, p, j)[1]:
        printf("Betting %s\n", convert(bet, string)):
        if GenerateFloat64() <= p then:
            cap := cap + bet:
            printf("You win! You have %d.\n", cap):
        else:
            cap := cap - bet:
            printf("You lose! You have %d.\n", cap):
        fi:
        if cap = N then:
            printf("You pay off the loan shark!"):
            return true:
        elif cap = 0 then:
            printf("You have run out of money, but at least you have %d days to live.", j):
            return false:
        fi:
    od:
    printf("As the time expires, the loan shark removes your head and places it on the ground."):
    return false:
end:

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

OptimalGamblingTerse := proc(N, i, p, T) local cap, j, bet, L, prob:
    cap := i:
    L := []:
    for j from T to 1 by -1 do:
        bet, prob := op(BestSDD(N, cap, p, j)):
        L := [op(L), prob]:
        if GenerateFloat64() <= p then:
            cap := cap + bet:
        else:
            cap := cap - bet:
        fi:
        if cap = N then:
            return true, L:
        elif cap = 0 then:
            # out of cash
            return false, L:
        fi:
    od:
    # out of time
    return false, L:
end:

#############
# Problem 2 #
#############
GetMin := proc(S, f:=(x->x)):
    return GetMax(S, x->-f(x)):
end:

ManyOptimalGambling := proc(N, i, p, T, K) local result, probs, c, L, j:
    c := 0:
    L := []:
    for j from 1 to K do:
        result, probs := OptimalGamblingTerse(N, i, p, T):
        if result then:
            c := c + 1:
            L := [op(L), GetMin(probs)[1]]:
        fi:
    od:
    return c/K, GetMin(L)[1]:
end:

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

# Win probability was 566/1000 which is within 0.001 of the true
# probability.

# Closest call was 0.073.





# More from my c24.txt:

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 4 #
#############

# The formulation in the homework makes no sense. Kelly is supposed
# to optimize the growth rate (i.e. the expected value of the capital
# after n steps is as large as possible, for every n). That is not the same as
# paying a unit cost per round.

# Nevertheless...

BestKellyFactor := proc(N, i, p, cost, res):
    return GetMax([seq(res..1, res)],
                  y->(VDS(p, N, x->1)[i+1]-VDS(p, N, x->ceil(y*x))[i+1])*cost - (VWS(p, N, x->1)[i+1]-VWS(p, N, x->ceil(y*x))[i+1])*N)[1]:
end:

#############
# Problem 5 #
#############

WhatKellyCharges := proc(p, N, i):
    if p <= 1/2 then:
        return 0:
    else:
        return (VWS(p, N, x->1)[i+1] - VWS(p, N, x->ceil((2*p-1)*x))[i+1])*N/(VDS(p, N, x->1)[i+1] - VDS(p, N, x->ceil((2*p-1)*x))[i+1]):
    fi:
end:



# WhatKellyCharges(11/20, 20, 10) = 0.0376
# WhatKellyCharges(11/20, 40, 21) = 0.0274
# WhatKellyCharges(11/20, 400, 150) = crashed maple

# I maintain that these values have no real significance.