game <- function(strategy, wealth, betsLeft) {
    if (betsLeft>0) {
    bet <- strategy(wealth, betsLeft)
    wealth <- wealth - bet
    flip <- rbinom(1,1,p=0.6)
    winnings <- 2*bet*flip
    wealth <- min(wealth+winnings, 250)
    return(game(strategy, wealth, betsLeft-1)) } else { return(wealth); } }
simulateGame <- function(s, w=25, b=300, iters=100) { mean(replicate(iters, game(s, w, b))) }

kelly <- function(w, b) { 0.20 * w }
smarterKelly <- function(w, b) { if(w==250) {0} else { (2*0.6-1) * w } }
random <- function(w, b) { sample(seq(0, w, by=0.1), 1) }

evaluate <- function(w,b) { simulateGame(smarterKelly, w, b) }
mevaluate <- memoise(evaluate)

mevaluate(25, 10)
# [1] 35.30544906
mevaluate(25, 100)
# [1] 159.385275
mevaluate(25, 300)
# [1] 231.4763619

f <- function(x, w, b, maxDepth) { 0.6*mV(min(w+x,250), b-1, maxDepth) + 0.4*mV(w-x, b-1, maxDepth)}
mf <- memoise(f)
V <- function(w,b, maxDepth) {
    if (b<=(b-maxDepth)) { mevaluate(w,b) } else {
    returns <- if (b>0 && w>0 && w<250) { sapply(seq(0, w, by=0.1), function(x) { mf(x,w,b, maxDepth) }) } else { w }
    max(returns) }}
mV <- memoise(V)

mV(25, 300, 100); mV(25, 300, 50)

(!!25) $ (!!300) $ (`iterate` [0..250]) $ \nextutilities -> (0:) $ (++[250]) $
    tail $ init $ map maximum $ (zipWith . zipWith) (\up down -> 0.6 * up + 0.4 * down)
     (tail $ tails nextutilities) (map reverse $ inits nextutilities)
-- 246.2080494949234
-- it :: (Enum a, Fractional a, Ord a) => a
-- (7.28 secs, 9,025,766,704 bytes)

:module + Control.Applicative
(!!25) $ (!!300) $ (`iterate` [0..250]) $ map maximum . liftA2 ((zipWith . zipWith)
    (\up down -> 0.6 * up + 0.4 * down)) tails (map reverse . tail . inits)
- 246.2080494949234
-- it :: (Enum a, Ord a, Fractional a) => a
-- (8.35 secs, 6,904,047,568 bytes)

import System.Environment (getArgs)
import Data.Function (iterate)
import Data.Function.ArrayMemoize (arrayMemo)

type Wealth = Int
type EV = Double

cap :: Wealth
cap = 250

value :: [Wealth -> EV]
value = iterate f fromIntegral where
  f next = arrayMemo (0, cap) $ \x -> maximum [go w | w <- [0..min (cap - x) x]] where
    go w = 0.6 * next (min cap (x+w)) + 0.4 * next (x-w)

main :: IO ()
main = do [x, b] <- map read <$> getArgs
          print $ value !! b $ x

import System.Environment (getArgs)
import Data.Function (fix)
import Data.Function.Memoize (memoFix2)
import Data.Vector (Vector)
import qualified Data.Vector as V (generate, (!))
import qualified Data.MemoUgly as U (memo)

type Wealth = Int
type Bets = Int
type EV = Double

cap :: Wealth
cap = 250

value :: (Wealth -> Bets -> EV) -> Wealth -> Bets -> EV
value next 0 b = 0
value next w 0 = fromIntegral w
value next w b = maximum [go x | x <- [0..min (cap - w) w]]
  where
    go x = 0.6 * next (min cap (w+x)) (b-1) + 0.4 * next (w-x) (b-1)

-- There are 4 possible ways aside from 'array-memoize':
-- 1. Direct recursion.
direct :: Wealth -> Bets -> EV
direct = fix value

-- 2. Memoized recursion.
memo :: Wealth -> Bets -> EV
memo w b = cached_value w b
  where
    cached_value w b = table V.! b V.! w
    table :: Vector (Vector EV)
    table = V.generate (b + 1)
            (\b -> V.generate (cap + 1)
              (\w -> value cached_value w b))

-- 3. Memoized recursion using 'memoize' library; slower.
memo2 :: Wealth -> Bets -> EV
memo2 = memoFix2 value

-- 4. Memoize using 'uglymemo' library; also slower but global
memo3 :: Wealth -> Bets -> EV
memo3 = U.memo . direct

main :: IO ()
main = do [w, b] <- map read <$> getArgs
          print (memo w b)

memo3 :: [Wealth] -> [Bets] -> [EV]
memo3 ws bs = zipWith cached_value ws bs
  where
    cached_value w b = table V.! b V.! w
    table :: Vector (Vector EV)
    table = V.generate (maximum bs + 1)
            (\b -> V.generate (cap + 1)
              (\w -> value cached_value w b))

λ> map round (memo3 (repeat 25) [1..300])
-- [30,36,43,45,48,54,55,58,62,63,67,68,71,74,75,79,80,82,85,86,90,90,94,94,97,99,100,
-- 103,104,107,108,110,111,113,115,116,119,120,122,123,125,126,128,130,131,133,134,
-- 136,136,139,139,141,142,144,145,146,148,149,151,151,153,154,155,156,158,159,160,
-- 161,162,164,164,166,166,168,169,170,171,172,173,174,175,176,177,177,179,179,181,
-- 181,182,183,184,185,186,186,187,188,189,190,190,191,192,193,193,194,195,196,196,
-- 197,198,198,199,200,200,201,202,202,203,204,204,205,205,206,207,207,208,208,209,
-- 209,210,210,211,212,212,213,213,214,214,214,215,215,216,216,217,217,218,218,219,
-- 219,219,220,220,221,221,221,222,222,222,223,223,224,224,224,225,225,225,226,226,
-- 226,227,227,227,228,228,228,228,229,229,229,230,230,230,230,231,231,231,231,232,
-- 232,232,232,233,233,233,233,234,234,234,234,234,235,235,235,235,236,236,236,236,
-- 236,236,237,237,237,237,237,238,238,238,238,238,238,239,239,239,239,239,239,239,
-- 240,240,240,240,240,240,240,241,241,241,241,241,241,241,241,242,242,242,242,242,
-- 242,242,242,242,243,243,243,243,243,243,243,243,243,243,244,244,244,244,244,244,
-- 244,244,244,244,244,244,245,245,245,245,245,245,245,245,245,245,245,245,245,245,
-- 246,246,246,246,246,246,246,246,246,246,246,246,246]

#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <math.h>

#define ROUND_LIMIT 300
#define MONEY_LIMIT 25000
#define MONEY_START 2500
#define MONEY_FACTOR 100.0f

// #define DATATYPE float
// #define F(X) X ## f
#define DATATYPE double
#define F(X) X

#define THREADS 8

// MEMOIZED METHOD

typedef struct {
  bool known;
  DATATYPE ev;
} CacheEntry;

static DATATYPE decide(CacheEntry * const table, int money, int round) {
  if (round < 0) return money;
  if (money == 0 || money == MONEY_LIMIT) return money;
  int index = round * (MONEY_LIMIT + 1) + money;
  // int index = money * ROUND_LIMIT + round;

  if (table[index].known) return table[index].ev;

  DATATYPE best_bet_score = -1;
  for (int bet = 0; bet <= money; bet++) {
    DATATYPE winres = decide(table, (money + bet > MONEY_LIMIT) ? MONEY_LIMIT : (money + bet), round - 1);
    DATATYPE loseres = decide(table, money - bet, round - 1);
    DATATYPE combined = F(0.6) * winres + F(0.4) * loseres;
    best_bet_score = F(fmax)(best_bet_score, combined);
  }
  table[index] = (CacheEntry) { true, best_bet_score };
  return best_bet_score;
}

void method1() {
  CacheEntry *table = calloc(sizeof(CacheEntry), ROUND_LIMIT * (MONEY_LIMIT + 1));

  printf("running.\n");
  int best_bet = 0;
  DATATYPE best_bet_score = -1;
  for (int bet = 0; bet <= MONEY_START; bet++) {
    DATATYPE winres = decide(table, MONEY_START + bet, ROUND_LIMIT - 1);
    DATATYPE loseres = decide(table, MONEY_START - bet, ROUND_LIMIT - 1);
    DATATYPE combined = F(0.6) * winres + F(0.4) * loseres;
    if (combined > best_bet_score) {
      best_bet = bet;
      best_bet_score = combined;
    }
  }

  printf("first round: bet %f for expected total reward of %f\n", best_bet / MONEY_FACTOR, best_bet_score / MONEY_FACTOR);
  int count_known = 0;
  for (int i = 0; i < ROUND_LIMIT * (MONEY_LIMIT + 1); i++) {
    if (table[i].known) count_known ++;
  }
  printf("known: %i of %i\n", count_known, ROUND_LIMIT * (MONEY_LIMIT + 1));
}

// BOTTOM-UP METHOD

// given pass n in "from", compute pass n-1 in "to"
static void propagate(DATATYPE *from, DATATYPE *to) {
  // for each money state...
#pragma omp parallel for num_threads(THREADS)
  for (int a = 0; a < THREADS; a++) {
    // distribute workload so that inner loop passes are approximately equal
    // f(x) = 2x  F(x) = x^2   = a/THREADS x = sqrt(a / THREADS)
    int low = (int) (MONEY_LIMIT * sqrtf((float)a / THREADS));
    int high = (int) (MONEY_LIMIT * sqrtf((a + 1.0f) / THREADS));
    if (a == THREADS - 1) high = MONEY_LIMIT + 1; // precise upper border
    for (int m = low; m < high; m++) {
      DATATYPE max_score = 0.0;
      // determine the bet that maximizes our expected earnings
      // to enable vectorization we break the loop into two so that each half only gets one branch of the conditional
      // DATATYPE winval = from[(m + b > MONEY_LIMIT) ? MONEY_LIMIT : (m + b)];
      int high_inner = (m + m > MONEY_LIMIT) ? (MONEY_LIMIT - m) : m;
      for (int b = 0; b <= high_inner; b++) {
        DATATYPE winval = from[m + b];
        DATATYPE loseval = from[m - b];
        DATATYPE combined = F(0.6) * winval + F(0.4) * loseval;
        max_score = F(fmax)(max_score, combined);
      }
      for (int b = high_inner + 1; b <= m; b++) {
        DATATYPE winval = from[MONEY_LIMIT];
        DATATYPE loseval = from[m - b];
        DATATYPE combined = F(0.6) * winval + F(0.4) * loseval;
        max_score = F(fmax)(max_score, combined);
      }
      to[m] = max_score;
    }
  }
}

void method2() {
  DATATYPE *buf = malloc(sizeof(DATATYPE) * (MONEY_LIMIT + 1));
  DATATYPE *buf_to = malloc(sizeof(DATATYPE) * (MONEY_LIMIT + 1));
  // set up base case: making no bet, we have money of i
  for (int i = 0; i <= MONEY_LIMIT; i++) buf[i] = i;
  for (int i = ROUND_LIMIT - 1; i >= 0; i--) {
    propagate(buf, buf_to);
    DATATYPE *temp = buf;
    buf = buf_to;
    buf_to = temp;
  }
  int best_bet = 0;
  DATATYPE best_bet_score = -1;
  for (int b = 0; b <= MONEY_START; b++) {
    DATATYPE winval = buf[MONEY_START + b];
    DATATYPE loseval = buf[MONEY_START - b];
    DATATYPE combined = 0.6 * winval + 0.4 * loseval;
    if (combined > best_bet_score) {
      best_bet_score = combined;
      best_bet = b;
    }
  }
  printf("first round: bet %f for expected total reward of %f\n", best_bet / MONEY_FACTOR, best_bet_score / MONEY_FACTOR);
}

int main() {
  method2();
  return 0;
}

for C in `ls coingame-*.c`; do
    gcc -Wall -Ofast -fwhole-program -march=native -fopenmp $C -o coingame -lm &&
    echo "$C" &&
    (for i in {1..10}; do
        /usr/bin/time -f '%e' ./coingame; done
    );
done

#include <iostream>
#include <vector>
#include <map>
#include <set>
#include <bitset>
#include <list>
#include <stack>
#include <queue>
#include <deque>
#include <string>
#include <sstream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <future>
#include <ctime>

using namespace std;

const int MAX_ROUND = 300;
const int MAX_WEALTH = 25000;
const int BETS_DELTA = 1;

double memo[MAX_WEALTH + 1][2];

inline double vWin(int wealth, int bet, int roundIdx) {
    if(wealth + bet < MAX_WEALTH) {
        return memo[wealth + bet][!roundIdx];
    }
    else {
        return MAX_WEALTH;
    }
}

inline double vLose(int wealth, int bet, int roundIdx) {
    if (wealth - bet == 0) {
        return 0;
    }
    else {
        return memo[wealth - bet][!roundIdx];
    }
}

inline double v(int wealth, int bet, int roundIdx) {
    return .6 * vWin(wealth, bet, roundIdx) + .4 * vLose(wealth, bet, roundIdx);
}

template<typename RAIter>
void calc_round(RAIter beg, RAIter end, int roundIdx){
    int len = distance(beg, end);

    if(len < 1000) {
        for(RAIter p = beg; p != end; ++p) {
            int wealth =  distance(memo, p);
            (*p)[roundIdx] = v(wealth, wealth, roundIdx);

            for (int bet = 0; bet < wealth; bet += BETS_DELTA) {
                memo[wealth][roundIdx] = max(memo[wealth][roundIdx], v(wealth, bet, roundIdx));
            }
        }
    }
    else {
        RAIter mid = beg + len/2;
        future<void> handle = async(launch::async, calc_round<RAIter>, mid, end, roundIdx);
        calc_round(beg, mid, roundIdx);
    }
}

double calc_table() {
    bool roundIdx = 0;

    for(int i = 0; i <= MAX_WEALTH; ++i) {
        memo[i][!roundIdx] = i;
    }

    for(int round = MAX_ROUND - 1; round >= 0; --round) {
        calc_round(memo, memo + MAX_WEALTH, roundIdx);
        roundIdx ^= 1;
    }

    return memo[2500][!roundIdx] / 100.;
}

int main() {
    ios_base::sync_with_stdio(0);
    cin.tie(0);
    clock_t begin = clock();
    double res = calc_table();
    clock_t end = clock();
    cout << "result: " << res << "(elapsed: " << (double(end - begin) / CLOCKS_PER_SEC) << ")" << endl;
}

V <- function(w, b, m=250) {
    b <- round(b) # rounds left must be integer
    j <- qbinom(w/m, b, 0.5)
    1.2^b * 1.5^-j * (w + (m/2 * sum(1.5^(j:0) * pbinom(0:j-1,b,1/2)))) }

from scipy.stats import binom
import numpy as np
from repoze.lru import lru_cache

def V(w, b, m=250):
    if w>=250:
      return 250
    if w<=0:
      return 0
    if b==0:
      return w
    else:
      try:
        j = binom.ppf(float(w)/float(m), b, 0.5)
        return 1.2**b * 1.5**-j * (w + m/2 *
            sum(np.multiply(binom.cdf(map(lambda(x2):x2-1, range(0,int(j+1))),b,0.5),
                map(lambda(x): 1.5**x, list(reversed(range(0, int(j+1))))))))
      except ValueError:
        print "Error:", (w,b,m)

@lru_cache(None)
def VPplan(w, b):
    # optimization: short-circuit
    if w<=0 or w>=250:
      return 0
    else:
      if b==0:
        return w
      else:
          possibleBets = map(lambda(pb): float(pb)/100.0, range(0*100,int((w*100)+1),1))
          returns = map(lambda(pb): 0.6*V(w+pb, b-1) + 0.4*V(w-pb,b-1), possibleBets)
          return float(returns.index(max(returns)))/100.0

import gym
env = gym.make('KellyCoinflip-v0')

## play 500 games and calculate mean reward:
rewards = []
for n in range(0,500):
    done = False
    reward = 0
    while not done:
        w = env._get_obs()[0][0]
        b = env._get_obs()[1]
        bet = VPplan(w, b)
        results = env._step(bet*100)
        print n, w, b, bet, "results:", results
        reward = reward+results[1]
        done = results[2]
    rewards.append(reward)
    env._reset()

print sum(rewards)/len(rewards)
# 247.525564356

# return all actions' values instead of the max action:
VPplans <- function(w, b) {
    if (b==0) { return (0); } else {
    returns <- sapply(seq(0, w), function(wp) { mf(wp, w, b); })
      return(returns); } }
plans <- data.frame()
n <- 47
for (i in 1:n) {
  # omit the uninteresting cases of w=0/250 and b=0/300:
  w <- round(runif(n=1, 0+1, 250-1))
  b <- round(runif(n=1, 0+1, 300-1))
  plan <- VPplans(w,b)
  plans <- rbind(plans, data.frame(W=as.factor(w), B=as.factor(b), Bet=1:(length(plan)), Value=plan))
}
library(ggplot2)
p <- qplot(x=Bet, y=Value, data=plans)
p + facet_wrap(~W) + geom_line(aes(x=Bet, y=Value, color=W), size = 1) + theme(legend.position = "none")

n <- 50000
w <- runif(n, min=0, max=250)
b <- sample(0:250, n, replace=TRUE)
v <- unlist(Map(V, w, b))
df <- data.frame(Wealth=w, Rounds=b, Value=v)

library(akima)
library(plot3D)
dfM <- interp(df$Wealth, df$Rounds, df$Value, duplicate="strip")
par(mfrow=c(1,2), mar=0)
persp(dfM$x, dfM$y, dfM$z, xlab="Wealth", ylab="Rounds left", zlab="Value", ticktype="detailed")
scatter3D(df$Wealth, df$Rounds, df$Value, theta=-25, xlab="Wealth", ylab="Rounds left", zlab="Value", ticktype="detailed")

library(randomForest)
r <- randomForest(Value ~ Wealth + Rounds, ntree=1000, data=df); r
# ...
#           Mean of squared residuals: 1.175970325
#                     % Var explained: 99.94
sqrt((predict(r) - df$Value)^2)
# [1] 0.407913389

## very slow
library(neuralnet)
n <- neuralnet(Value ~ Wealth + Rounds, data=df, hidden=c(20,20))
# prediction(n)

## GPU-accelerated:
library(mxnet)
mx.set.seed(0)
## 2×30 ReLu FC NN for linear regression with RMSE loss:
data <- mx.symbol.Variable("data")
perLayerNN <- 30
fc1     <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=perLayerNN)
act1    <- mx.symbol.Activation(fc1, name="relu1", act_type="relu")
fc2     <- mx.symbol.FullyConnected(act1, name="fc2", num_hidden=perLayerNN)
act2    <- mx.symbol.Activation(fc2, name="relu2", act_type="relu")
fc3     <- mx.symbol.FullyConnected(act2, name="fc2", num_hidden=perLayerNN)
act3    <- mx.symbol.Activation(fc3, name="relu2", act_type="relu")
finalFc <- mx.symbol.FullyConnected(act3, num_hidden=1)
lro <- mx.symbol.LinearRegressionOutput(finalFc)
model <- mx.model.FeedForward.create(lro, X=as.matrix(subset(df,select=c("Wealth","Rounds"))), y=df$Value,
    ctx=mx.gpu(),     num.round=50000,
    learning.rate=5e-6, eval.metric=mx.metric.rmse, array.layout="rowmajor")