heads-up

#!/usr/bin/env python

from __future__ import division
import os
import re
import sys
import optparse
from numpy import *
from nash import *
from util import *

usage = "usage: %prog [options...] [args...]"
parser = optparse.OptionParser(usage)
parser.add_option('-b','--bet',type=str,default=0,help='set a fixed bet level')
parser.add_option('--bets',type=str,default=[],help='explicitly set the possible bets')
parser.add_option('--sb',type=float,default=1,help='small blind')
parser.add_option('--bb',type=float,default=2,help='big blind')
parser.add_option('--plot',action='store_true',help="plot Alice's equity for bets in [0,bet]")
parser.add_option('--max',action='store_true',help="maximize Alice's equity as a function of bet in [0,bet]")
parser.add_option('--donkey',type=int,default=0,help="compute optimum strategy against an always all-in player, with a given total chip count")
parser.add_option('--random',action='store_true',help="print win/lose/tie probabilities of each hand vs. a random hand")
parser.add_option('--variable',type=int,default=0,help="compute the Nash equilibrium for a game where Alice has a choice of bets")
parser.add_option('--inexact',action='store_true',help='avoid use of exact rational arithmetic')
parser.add_option('-d','--deterministic',action='store_true',help="switch to deterministic outcomes")
parser.add_option('-n',type=int,default=10,help="number of samples to plot")
parser.add_option('-v','--verbose',action='store_true',help="print or plot extra information")
parser.add_option('-p','--profile',action='store_true',help="skip plotting for profiling purposes")
parser.add_option('-t','--type',type=str,help="type of plot to show")
options,args = parser.parse_args()
if len(args): parser.error('zero arguments expected')

set_printoptions(linewidth=150)

if not options.inexact:
  from rational import *
  R = rational
  def veryclose(a,b):
    return all(a==b)
else:
  R = float
  veryclose = allclose

cards = '23456789TJQKA'
nonpairs = [cards[i]+cards[j] for i in xrange(len(cards)) for j in xrange(i)]
hands = [c+c for c in cards]+[h+'s' for h in nonpairs]+[h+'o' for h in nonpairs]
assert len(hands)==13+13*12
hand_id = dict((h,i) for (i,h) in enumerate(hands))

def die(s):
    print>>sys.stderr, s
    sys.exit(1)

def ahash(*args):
    return hash(''.join(asarray(a).tostring() for a in args))

# Load exact win probabilities.  Important: win[a,b]+win[b,a]<1 due to ties, unlike the approximate version above
def parse_exact_matchups():
    file = 'exact.npz'
    if not os.path.exists(file):
        matchup_pattern = re.compile(r'^(\w+)\s+vs\.\s+(\w+):$')
        entry_pattern = re.compile(r'^\s+(Alice|Bob|Tie):\s+(\d+/\d+) = ([\d.]+)$')
        lines = open('exact.txt').readlines()
        assert len(lines)%4==0
        win = zeros((len(hands),)*2,dtype=R)
        for i in xrange(len(lines)//4):
            m = matchup_pattern.match(lines[4*i])
            if not m: die('Weird matchup line: %s'%lines[4*i][:-1])
            h0,h1 = hand_id[m.group(1)],hand_id[m.group(2)]
            def entry(j,name):
                m = entry_pattern.match(lines[4*i+j+1])
                if not m: die('Weird entry line: %s'%lines[4*i+j+1][:-1])
                assert m.group(1)==name
                if options.inexact:
                  a,b = map(float64,m.group(2).split('/'))
                  return a/b
                return rational(m.group(2))
            w,l,t= entry(0,'Alice'),entry(1,'Bob'),entry(2,'Tie')
            assert veryclose(w+l+t,1)
            win[h0,h1] = w
            win[h1,h0] = l
        assert all(win>0)
        if options.inexact:
          savez_compressed(file,win=win)
        else:
          savez_compressed(file,n=numerator(win),d=denominator(win))
        return win
    else:
        data = load(file)
        if options.inexact:
          win = data['win']
        else:
          win = rationals(data['n'])/data['d']
        return win
exact_win = parse_exact_matchups()
assert ahash(exact_win)==(8198859743626163791 if options.inexact else 7392642596451405510)

# Build a map from pairs of cards (represented as 4*id+suit) to hand ids
def compute_cards_to_hand():
    suits = xrange(4)
    cards_to_hand = -ones((4*len(cards),)*2,dtype=int)
    for c0 in xrange(len(cards)):
        h = hand_id[cards[c0]*2]
        for s0 in suits:
            for s1 in suits:
                if s0!=s1:
                    cards_to_hand[4*c0+s0,4*c0+s1] = h
        for c1 in xrange(c0):
            h = cards[c0]+cards[c1]
            hs = hand_id[h+'s']
            ho = hand_id[h+'o']
            for s0 in suits:
                cs0 = 4*c0+s0
                for s1 in suits:
                    cs1 = 4*c1+s1
                    cards_to_hand[cs0,cs1] = cards_to_hand[cs1,cs0] = hs if s0==s1 else ho
    return cards_to_hand
cards_to_hand = compute_cards_to_hand()
assert sum(cards_to_hand<0)==52

def hand_probabilities():
    # Count how many times h occurs in cards_to_hand for each h
    def prob(x):
        prob = sum(x.reshape(1,-1)==arange(len(hands)).reshape(-1,1),axis=1)
        return prob.astype(R)/sum(prob)
    hand_prob = prob(cards_to_hand)

    i = xrange(len(cards_to_hand))
    hand_to_cards = dict((cards_to_hand[c,d],(c,d)) for c in i for d in i)

    # Compute conditional probabilities of each hand given each other hand
    cond_hand_prob = zeros((len(hands),)*2,R)
    for h in xrange(len(hands)):
        cs = hand_to_cards[h]
        remaining = array([i for i in xrange(13*4) if i not in cs])
        cond_hand_prob[h] = prob(cards_to_hand[remaining.reshape(-1,1),remaining.reshape(1,-1)])
    return hand_prob,cond_hand_prob
exact_hand_prob,exact_cond_hand_prob = hand_probabilities()
exact_hand_hand_prob = exact_hand_prob.reshape(-1,1)*exact_cond_hand_prob
assert veryclose(sum(exact_cond_hand_prob,axis=1),1)
assert veryclose(sum(exact_hand_hand_prob),1)
assert veryclose(exact_hand_hand_prob,exact_hand_hand_prob.T)

# Check all interesting conditional probabilities
assert veryclose(exact_hand_prob[hand_id['AA']],R(1)/13*3/51)
assert veryclose(exact_hand_prob[hand_id['AKs']],R(2)/13/51)
assert veryclose(exact_hand_prob[hand_id['AKo']],R(2)/13*3/51)
assert veryclose(exact_cond_hand_prob[hand_id['AA'],hand_id['AA']],R(2)/50*1/49)
assert veryclose(exact_cond_hand_prob[hand_id['AA'],hand_id['KK']],R(4)/50*3/49)
assert veryclose(exact_cond_hand_prob[hand_id['AA'],hand_id['AKs']],R(2)*2/50/49)
assert veryclose(exact_cond_hand_prob[hand_id['AA'],hand_id['AKo']],R(2)*2/50*3/49)
assert veryclose(exact_cond_hand_prob[hand_id['AA'],hand_id['KQs']],R(2)*4/50/49)
assert veryclose(exact_cond_hand_prob[hand_id['AA'],hand_id['KQo']],R(2)*4/50*3/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKs'],hand_id['AKs']],R(2)*3/50/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKs'],hand_id['AKo']],R(2)*3/50*2/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKs'],hand_id['AQs']],R(2)*3/50/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKs'],hand_id['AQo']],R(2)*3/50*3/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKs'],hand_id['QJs']],R(2)*4/50/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKs'],hand_id['QJo']],R(2)*4/50*3/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKo'],hand_id['AKo']],R(2)/50*3/49+R(2)*2/50*2/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKo'],hand_id['AQo']],R(2)*3/50*3/49)
assert veryclose(exact_cond_hand_prob[hand_id['AKo'],hand_id['QJo']],R(2)*4/50*3/49)

# Switch to deterministic outcomes
if options.deterministic:
    outcome = dot(exact_win,exact_hand_prob)
    exact_win[...] = outcome.reshape(-1,1)>outcome.reshape(1,-1)
    i = arange(len(outcome))
    assert all(0<=exact_win) and all(exact_win<=1)
    rank = sorted(i,key=lambda i:outcome[i])
    print 'hand rank = %s'%(' '.join(hands[i] for i in rank))

# Approximate
approx_win = exact_win.astype(float)
approx_hand_prob = exact_hand_prob.astype(float)
approx_hand_hand_prob = exact_hand_hand_prob.astype(float)
approx_cond_hand_prob = exact_cond_hand_prob.astype(float)

def if_(c,a,b):
    return c*a+(1-c)*b # Use multiplication in order to work for all probabilities in [0,1]

def check_payoff(bet,alice,bob,exact=True):
    if exact:
        hand_hand_prob,win = exact_hand_hand_prob,exact_win
    else:
        bet = float(bet)
        hand_hand_prob,win = approx_hand_hand_prob,approx_win
    return sum(hand_hand_prob*alice.reshape(-1,1)*if_(bob.reshape(1,-1),(1+bet)*(win-win.T),1))

def alice_payoff(bet,bob,exact=True):
    "Alice's expected return from calling with each possible hand"
    if exact:
        cond_hand_prob,win = exact_cond_hand_prob,exact_win
    else:
        bet = float(bet)
        cond_hand_prob,win = approx_cond_hand_prob,approx_win
    return sum(cond_hand_prob*if_(bob.reshape(1,-1),(1+bet)*(win-win.T),1), axis=1)

def bob_payoff(bet,alice,exact=True):
    "Bob's expected return from calling with each possible hand"
    if exact:
        hand_hand_prob,win = exact_hand_hand_prob,exact_win
    else:
        bet = float(bet)
        hand_hand_prob,win = approx_hand_hand_prob,approx_win
    hand_hand_call_prob = alice.reshape(1,-1)*hand_hand_prob
    inv_cond_hand_prob = hand_hand_call_prob / hand_hand_call_prob.sum(axis=1).reshape(-1,1)
    return sum(inv_cond_hand_prob*(1+(1+bet)*(win-win.T)), axis=1)

def always_call():
    alice = ones(len(hands),dtype=int)
    bob = ones(len(hands),dtype=int)
    return alice,bob

# Converge to heads up Nash equilibrium from an initial strategy guess
def poker_nash_equilibrium(bet,(alice,bob)=always_call(),n=1000,verbose=1):
    alice = alice.astype(bool)
    bob = bob.astype(bool)
    bounds = [-inf,inf]
    # Iterate until we detect looping
    hashes = set()
    looping = False
    strategies = []
    while 1:
        done = 1
        if verbose:
            print 'equity range = %g %g'%(bounds[0],bounds[1])
        # Make at most one change to Alice's strategy
        payoff = alice_payoff(bet,bob,exact=0)
        bounds[1] = min(bounds[1],dot(approx_hand_prob,maximum(0,payoff))) # Bob can hold Alice to below this
        error = payoff*((payoff>0)*(1-alice) - (payoff<0)*alice)
        h = argmax(error)
        if error[h]:
            alice = alice.copy()
            alice[h] = 1-alice[h]
            if verbose:
                print 'Alice should%s call with %s'%(('' if alice[h] else "n't"),hands[h])
            done = 0
        # Make at most one change to Bob's strategy
        payoff = bob_payoff(bet,alice,exact=0)
        error = payoff*((payoff>0)*(1-bob) - (payoff<0)*bob)
        bounds[0] = max(bounds[0],check_payoff(bet,alice,payoff>0,exact=0)) # Alice can hold Bob to above this
        h = argmax(error)
        if error[h]:
            bob = bob.copy()
            bob[h] = 1-bob[h]
            if verbose:
                print 'Bob should%s call with %s'%(('' if bob[h] else "n't"),hands[h])
            done = 0
        # Are we done?
        if done:
            strategies = [(alice,bob)]
            break
        ah = ahash(alice,bob)
        if verbose:
            print 'hash = %s'%ah
        if ah in hashes:
            if looping:
                break
            else:
                hashes = set()
            looping = True
        if looping:
            strategies.append((alice,bob))
        hashes.add(ah)
    # Print strategies
    print '\nfinal equity range = %g %g'%(bounds[0],bounds[1])
    for k,name in enumerate(['Alice','Bob']):
        print '%s has %d %s:'%(name,len(strategies),'strategy' if len(strategies)==1 else 'strategies')
        for s in strategies:
            call = s[k]
            print '  %d hands: %s'%(sum(call),' '.join(hands[i] for i in nonzero(call)[0]))
    # Compute exact Nash equilibrium matrix
    alices = [s[0] for s in strategies]
    bobs = [s[1] for s in strategies]
    while 1:
        # Compute Nash equilibrium ranging over our current list of Alice's and Bob's strategies
        print 'computing payoffs...'
        payoff = empty((len(alices),len(bobs)),R)
        for a,alice in enumerate(alices):
            for b,bob in enumerate(bobs):
                payoff[a,b] = check_payoff(bet,alice,bob,exact=1)
        print 'computing payoffs...done'
        nash,alice,bob = zero_sum_nash_equilibrium(payoff)
        alice = sum(a*s for a,s in zip(alice,alices))
        bob = sum(b*s for b,s in zip(bob,bobs))
        # Determine Alice's and Bob's new optimal strategies, and add them to our lists if necessary
        done = 1
        opt_alice = alice_payoff(bet,bob)>0
        if ahash(opt_alice) not in map(ahash,alices):
            alices.append(opt_alice)
            done = 0
        opt_bob = bob_payoff(bet,alice)>0
        if ahash(opt_bob) not in map(ahash,bobs):
            bobs.append(opt_bob)
            done = 0
        if done:
            break
    # Verify that we've found the Nash equilibrium
    equity = check_payoff(bet,alice,bob)
    assert equity==nash
    assert bounds[0]-1e10<=equity<=bounds[1]+1e10
    if 0:
        print 'first error = %g'%(check_payoff(bet,alice,bob_payoff(bet,alice)>0)-equity)
        print 'second error = %g'%(check_payoff(bet,alice_payoff(bet,bob)>0,bob)-equity)
    assert equity==check_payoff(bet,alice,bob_payoff(bet,alice)>0)
    assert equity==check_payoff(bet,alice_payoff(bet,bob)>0,bob)
    # Print out nice results
    if verbose:
        def str_strategy(strat,payoff):
            order = sorted(arange(len(strat)),key=lambda i:-payoff[i])
            return ' '.join(('' if strat[i]==1 else '%s*'%strat[i])+hands[i] for i in order if strat[i])
        print
        print 'Alice: %d hands, %s'%(sum(alice>0),str_strategy(alice,alice_payoff(bet,bob)))
        print 'Bob: %d hands, %s'%(sum(bob>0),str_strategy(bob,bob_payoff(bet,alice)))
    return equity,alice,bob

def poker_nash_equilibrium_remember(bet,last=list(always_call())):
    equity,alice,bob = poker_nash_equilibrium(bet,last)
    last[:] = round_(alice),round_(bob)
    return equity

def donkey_strategy(total):
    """Game: Alice vs. Bob, alternating small blind of 1 and big blind of 2.  Bob always raises all in.
    What is Alice's optimal strategy given a total of T chips?  Let Alice's probability of winning be
        equity : {SB=0,BB=1} -> {0...T} -> [0,1]
    We have equity[s,0] = 0, equity[s,T] = T.  Let win[a,b] be the probability of hand a winning against hand b,
        win[a] = sum_b p[b|a]win[a,b]
    be the probability of winning against a random hand, and
        lose[a] = sum_b p[b|a]win[b,a]
    the probability of losing against a random hand.  If Alice folds, she loses 1+s, where s=0 for SB and
    1 for BB.  If Alice calls with a stack of t chips, then her outcomes are
        +min(t,T-t) with prob win[a]
        0           with prob 1-win[a]-lose[a] = tie[a]
        -min(t,T-t) with prob lose[a]
    for a total equity of 
        win[a]equity[1-s,min(2t,T)] + tie[a]equity[1-s,t] + lose[a]equity[1-s,max(0,2t-T)]
    Thus Alice's equity at equity[s,t] is
        equity[s,t] = sum_a p[a] max(equity[1-s,max(0,t-1-s)], win[a]equity[1-s,min(2t,T)] + tie[a]equity[1-s,t] + lose[a]equity[1-s,max(0,2t-T)])
    This is something we can iterate.
    """
    win = sum(approx_cond_hand_prob*approx_win,axis=1)
    lose = sum(approx_cond_hand_prob*approx_win.T,axis=1)
    tie = 1-win-lose
    equity = zeros((2,total+1))
    equity[:,total] = 1
    call = zeros((2,total+1,len(hands)),dtype=bool)
    s = arange(2).reshape(-1,1)
    t = arange(1,total).reshape(1,-1)
    def hand_reshape(x):
        return x.reshape((2,total-1,1))
    for iter in xrange(1000000):
        print 'strategizing...'
        equity_fold = equity[1-s,maximum(0,t-1-s)]
        equity_call = win*hand_reshape(equity[1-s,minimum(2*t,total)])+tie*hand_reshape(equity[1-s,t])+lose*hand_reshape(equity[1-s,maximum(0,2*t-total)])
        new_call = equity_call>hand_reshape(equity_fold)
        if all(call==new_call):
            print 'converged'
            break
        call = new_call

        print 'building...'
        import scipy.sparse
        import scipy.sparse.linalg
        ii,jj,vv = [],[],[]
        def add(i,j,v):
            i,j = i+0*j,j+0*i # broadcast
            if not ndim(v): v = v+0*i
            assert i.shape==j.shape==v.shape
            ii.append(ravel(i))
            jj.append(ravel(j))
            vv.append(ravel(v))
        add(s+2*t,s+2*t,-1)
        add(s+2*t,1-s+2*maximum(0,t-1-s),sum(approx_hand_prob.reshape(1,1,-1)*(1-call),axis=-1))
        prob_call = approx_hand_prob.reshape(1,1,-1)*call
        add(s+2*t,1-s+2*minimum(2*t,total),sum(prob_call*win,axis=-1))
        add(s+2*t,1-s+2*t,sum(prob_call*tie,axis=-1))
        add(s+2*t,1-s+2*maximum(0,2*t-total),sum(prob_call*lose,axis=-1))
        As = scipy.sparse.coo_matrix((hstack(vv),vstack([hstack(ii),hstack(jj)])),(2*total+2,2*total+2)).tocsr()
        new_equity = equity.copy()
        b = -As[2:-2,-2:].sum(axis=1)
        A = As[2:-2,2:-2]
        print 'solving...'
        new_equity[s,t] = scipy.sparse.linalg.spsolve(A,b).reshape(total-1,2).T

        assert all(-1e-7<new_equity)
        assert all(new_equity<1+1e-7)
        assert all(new_equity+1e-7>equity)
        equity = new_equity
        print "iteration %d, Alice's small blind equity with %d chips is %g"%(iter,total//2,equity[0,total//2])
    return equity,call

def limit_donkey_equity(k):
    """In the case of infinite stack size, Alice's optimal strategy is simply to wait for AA.  The resulting
    win probability function satisfies
        p[t] = 0, t <= 0
        p[t] = 1, t >= 1
        p[t] = a p[2t] + (1-a) p[2t-1], 0 < t < 1
    where a is the probability of winning with AA (splitting ties).  Important properties:
        1. p is continuous.
        2. If t is a dyadic rational, p[t] is a function of simpler dyadic rationals.
    This provides an easy, superfast way of computing p."""
    a = 893604787/1048786200
    p = zeros(1)
    for k in xrange(k):
        p = hstack([a*p,a+(1-a)*p])
    return hstack([p,1])

def variable_bet_nash_equilibrium(sb,bb,bets):
    """We consider a heads up game with Alice and Bob where

    1. Alice posts a small blind of sb.
    2. Bob posts a big blind of bb.
    3. Alice either folds or bets one of b in bets.
    4. Bob either folds or calls.

    The set of hands is H, the joint probability of hands x and y is pp(x,y), and the probability that x beats y is w(x,y).  Ignore ties so that w(x,y) + w(y,x) = 1.  I.e.,"""
    pp = approx_hand_hand_prob
    ps = pp.sum(axis=1)
    w = approx_win + (1-approx_win-approx_win.T)/2
    """Alice's strategy is a probability Pij of betting amount i given hand j.  Bob's strategy is a probability Qik of calling given Alice's bet i and Bob's hand k.
    We shift all outcomes by shift so that Alice's equity is uniformly positive:
    
        A = shift - sb sum_j pj (1-sum_i Pij) + sum_ijk pjk Pij (bb (1-Qik) + Qik (bb+bet[i])(wjk - wkj))
          = (shift-sb) + sb sum_ij pj Pij + sum_ijk pjk Pij (bb + Qik ((bb+bet[i])(wjk-wkj)-bb))
          = (shift-sb) + (sb+bb) sum_ij pj Pij + sum_ijk pjk Pij Qik ((bb+bet[i])(wjk-wkj)-bb)
          = mu + sum_ij gij Pij + sum_ijk Gijk Pij Qik

    where

        mu = shift-sb
        gij = (sb+bb) pj
        Gijk = pjk ((bb+bets[i])(wjk-wkj)-bb)
        shift = max(sb,bb)+max(sb,bb+bets.max())
    
    Let's construct mu, g and G:"""
    shift = max(sb,bb)+max(sb,bb+bets.max())
    mu = shift-sb
    g = (sb+bb)*ps.reshape(1,-1)*ones((len(bets),1))
    i = arange(len(bets)).reshape(-1,1,1)
    j = arange(len(hands)).reshape(1,-1,1)
    k = j.reshape(1,1,-1)
    G = pp[j,k]*((bb+bets[i])*(w[j,k]-w[k,j])-bb) # Flexible indexing is sweet
    """In order to write the equity as a standard quadratic form, let p = ravel(P), q = ravel(Q), and define a block matrix G by

        Gab = G(ij)(ikl) = Gijkl"""
    g = ravel(g)
    G = spdiag(G.reshape((len(bets),len(hands),len(hands))))
    """Since probabilities sum to one, we have sum_i Pij <= 1, = Qik <= 1.  The scale of P will have to vary, so add an extra
    scalar element u to p, and enforce sum_i Pij <= u instead.  In matrix form, these constraints are
        Bp <= 0, B(j,ij) = 1, B(j,-1) = -1
        q <= 1"""
    B = sparse.hstack([sparse.csr_matrix((ones(len(bets)*len(hands)),arange(len(bets)*len(hands)).reshape(len(bets),-1).T.ravel(),len(bets)*arange(len(hands)+1))),-ones((len(hands),1))])
    g = hstack([g,mu]) # Absorb mu into g
    G = sparse.vstack([G,zeros((1,len(bets)*len(hands)))])
    """Before scaling, the problem we need to solve is

        A = max_{p>=0, Bp<=0, u=1} min_{q>=0, q<=1} (mu + g'p + p' G q)

    Allowing the scale variable u to grow until g'p + p' G q >= 1, we get

        1/A = min u s.t. p>=0, Bp<=0 (inf_{q>=0, q<=1} p' G q)>=1-g'p

    Following Tomas Tinoco De Rubira's email, we will dualize the universally quantified inequality into an existentially qualified inequality.  The primal problem is

        min   (G'p)'q
        s.t.  q >= 0
        s.t.  q <= 1

    which dualizes to
    
        max   -1'z
        s.t.  z + G'p >= 0
              z >= 0

    Strong duality gives

        (inf_{q>=0, q<=1} p' G q) = (sup_{z>=0, G' p + z >= 0} -1'z)
        (1-g'p<=inf_{q>=0, q<=1} p^T G q) iff exists z>=0 s.t. g'p - 1'z >= 1, G' p + z >= 0

    and our problem becomes

        1/A = min u s.t. p>=0, z>=0, Bp<=0, g'p - 1'z >= 1, G'p + z >= 0
            = min u s.t. p>=0, z>=0, Bp<=0, -g'p + 1'z <= -1, -G'p - z <= 0

    Now write x = [p,z], and construct H = [-1,[-G.T,-1],[-g',1'],[B,0]], h = [0,0,-1,0]:"""
    pn,qn = len(bets)*len(hands)+1,len(bets)*len(hands)
    H = sparse.vstack([-speye(pn+qn),sparse.bmat([[-G.T,-speye(qn)],[sparse.coo_matrix(-g.reshape(1,-1)),ones((1,qn))],[B,None]])])
    h = hstack([zeros(pn+2*qn),-1,zeros(len(hands))])
    e = hstack([zeros(pn-1),1,zeros(qn)])
    """Our problem is now

        1/A = min e'x s.t. Hx <= h

    which is exactly the form expected by cvxopt.  Let's verify that the program is feasible:"""
    x = hstack([zeros(pn-1),1/mu,zeros(qn)])
    assert all(H*x-h<=1e-10)
    """Time to think for a while:"""
    solution = cvxopt_lp(e,H,h)
    assert solution['status']=='optimal'
    """Alice's optimal strategy and equity are"""
    x = solution['x']
    p,u,_ = asplit(x,pn-1,1,qn)
    P = p.reshape(len(bets),len(hands))/u
    equity = 1/u-shift
    print "Alice's equity = %g"%equity
    """What about Bob's Nash equilibrium strategy?  Can we extract it from the dual solution that we already have, rather than explicitly solving another dual LP?
    The dual to the program we gave to cvxopt is

        max  -h'y
        s.t. H'y + e = 0
             y >= 0
 
    where y = [y0,y1,y2,y3,y4].  Expanding H'y + e = 0 gives

        -y0 - G y2 - g y3 + B' y4 + e = 0
        -y1 - y2 + 1 y3 = 0

    or

        -G y2 - g y3 + B' y4 + e >= 0
        y2 <= 1 y3

    hinting that perhaps q = y2/y3.  We test this by constructing Alice's optimal strategy assuming Bob's strategy is q = y2/32:"""
    y = solution['z']
    _,q,qu,_ = asplit(y,pn+qn,qn,1,len(hands))
    Q = q.reshape(len(bets),len(hands))/qu
    return equity,P,Q

def check_variable_bet_nash_equilibrium(sb,bb,bets,equity,P,Q):
    assert all(P>=-1e-4)
    assert all(sum(P,axis=0)<1+1e-4)
    assert all(Q>=-1e-4)
    assert all(Q<=1+1e-4)
    """Recompute some variables from the above routine"""
    pp = approx_hand_hand_prob
    ps = pp.sum(axis=1)
    w = approx_win + (1-approx_win-approx_win.T)/2
    """Now we work out Bob's optimal strategy given Alice's fixed strategy P.  Given Alice's bet i and Bob's hand k, the conditional probability of Alice having hand j is
        ca[i,j] = pp[j,k] P[i,j] / sum_j pp[j,k] P[i,j]"""
    ca = pp.reshape(1,len(hands),len(hands))*P.reshape(len(bets),len(hands),1)
    ca /= ca.sum(axis=1).reshape(len(bets),1,len(hands))
    """If Bob folds, his equity is -bb.  If Bob calls, his equity is
        Be[i,k] = sum_j ca[i,j,k] (bb+bets[i]) (w[k,j]-w[j,k])"""
    Be = sum(ca*(bb+bets).reshape(-1,1,1)*(w.T-w),axis=1)
    equity2 = -sb*sum(pp*(1-P.sum(axis=0)).reshape(-1,1)) + sum(pp*P.reshape(len(bets),len(hands),1)*if_((Be<-bb).reshape(len(bets),1,len(hands)),bb,(bb+bets.reshape(-1,1,1))*(w-w.T)))
    Be = maximum(Be,-bb)
    print "Alice's equity if Bob switches = %g"%equity2
    print '  error = %g'%abs(equity-equity2)
    """Now we compute Alice's optimal strategy given Bob's fixed strategy Q."""
    cb = pp/sum(pp,axis=1).reshape(-1,1)
    Ae = empty((len(bets)+1,len(hands)))
    Ae[0] = -sb
    Ae[1:] = sum(cb*if_(Q.reshape(len(bets),1,len(hands)),(bb+bets.reshape(-1,1,1))*(w-w.T),bb),axis=-1)
    P2 = zeros((len(bets)+1,len(hands)))
    P2[argmax(Ae,axis=0),arange(len(hands))] = 1
    equity3 = -sb*sum(ps*P2[0]) + sum(pp*P2[1:].reshape(len(bets),len(hands),1)*if_(Q.reshape(len(bets),1,len(hands)),(bb+bets.reshape(-1,1,1))*(w-w.T),bb))
    assert allclose(dot(Ae.max(axis=0),approx_hand_prob),equity3)
    print "Alice's equity if Alice switches = %g"%equity3
    print '  error = %g'%abs(equity-equity3)
    print 'Equity spread = %g'%(equity3-equity2)
    """All done!"""
    return (equity2,equity3),Ae,Be

bet = options.bet
try:
    bet = rational(bet)
except TypeError:
    bet = float(bet) 

if options.plot:
    n = arange(options.n)
    bet = bet/n[-1]*n
    print 'bets = %s'%bet
    s = list(always_call())
    def nash(b):
        e,alice,bob = poker_nash_equilibrium(b,s,verbose=1)
        s[:] = round_(alice),round_(bob)
        return e
    equity = array(map(poker_nash_equilibrium_remember,bet))
    import pylab
    pylab.plot(bet,equity)
    pylab.xlabel('bet')
    pylab.ylabel("Alice's equity")
    pylab.title("Alice's equity as a function of fixed bet size")
    pylab.show()
elif options.max:
    import scipy.optimize
    bet,equity,_,evals = scipy.optimize.fminbound(lambda b:-poker_nash_equilibrium_remember(b),0,bet,full_output=1,xtol=1e-8)
    print "\nAlice's preferred bet = %g"%bet
    print "Alice's equity = %g"%-equity
    print "Nash equilibrium computations = %d"%evals
elif options.donkey>0:
    total = options.donkey
    equity,call = donkey_strategy(total)
    call_hash = ahash(call)
    expected = {200:8011577361726862008,2000:-5197011106140037419,20000:654298559373924011}.get(total,'unknown')
    print 'call hash: expected %s, got %d'%(expected,call_hash)
    if options.profile:
        sys.exit(0)
    call_prob = sum(approx_hand_prob.reshape(1,1,-1)*call,axis=-1)
    t = total//2
    print 'With %d chips in the small blind, Alice calls %g of the time, with equity %g'%(t,call_prob[0,t],equity[0,t])
    print 'With %d chips in the big blind, Alice calls %g of the time, with equity %g'%(t,call_prob[1,t],equity[1,t])
    import pylab
    t = arange(1,total)
    pylab.plot(t,equity[0,t],'g',label='SB equity')
    pylab.plot(t,equity[1,t],'b',label='BB equity')
    if options.verbose:
        win = sum(approx_cond_hand_prob*approx_win,axis=1)
        ordered_hands = argsort(win)
        hand_rank = empty_like(ordered_hands)
        hand_rank[ordered_hands] = arange(len(hands))
        if 0:
            print "Alice's calling hands:"
            for i in t:
                for s in 0,1:
                    print '  %3d %s, %s: %s'%(i,('chip' if i==1 else 'chips'),('SB','BB')[s],' '.join(hands[h] for h in sorted(nonzero(call[s,i-1])[0],key=lambda h:-win[h])))
        for s in 0,1:
            for t in xrange(1,total,3):
                min_rank = hand_rank[nonzero(call[s,t-1])[0]].min()
                pylab.text(t,call_prob[s,t-1],hands[ordered_hands[min_rank]],horizontalalignment='center',verticalalignment='center',size='xx-small',color='mr'[s])
        # For legend purposes
        pylab.plot([],[],'mo',label='SB worst call hand')
        pylab.plot([],[],'ro',label='BB worst call hand')
    else:
        pylab.plot(t,call_prob[0],'mo',label='SB call prob')
        pylab.plot(t,call_prob[1],'ro',label='BB call prob')
    pylab.legend(loc=('upper center' if total<=500 else 'center'))
    pylab.xlabel("Alice's stack size")
    pylab.ylabel("Equity or call probability")
    pylab.title('Optimal strategy vs. an always all-in player')
    pylab.show()
elif options.donkey<0:
    k = -options.donkey
    equity = limit_donkey_equity(k)
    t = arange(2**k+1)/2**k
    import pylab
    pylab.plot(t,equity,'b')
    pylab.xlabel("Alice's stack fraction")
    pylab.ylabel("Alice's equity")
    pylab.title('Optimal equity vs. an always all-in player with infinite stack sizes')
    pylab.show()
elif options.random:
    win = sum(exact_cond_hand_prob*exact_win,axis=1)
    lose = sum(exact_cond_hand_prob*exact_win.T,axis=1)
    tie = 1-win-lose
    ordered_hands = argsort(win)
    print 'Probabilities vs. a random hand:'
    for h in ordered_hands[::-1]:
        print '  %3s : win %s, lose %s, tie %s, win+tie/2 %s'%(hands[h],win[h],lose[h],tie[h],win[h]+tie[h]/2)
        #print '  %3s : win %g, lose %g, tie %g, win+tie/2 %g'%(hands[h],win[h],lose[h],tie[h],win[h]+tie[h]/2)
elif options.variable or options.bets:
    sb = options.sb
    bb = options.bb
    if options.bets:
        bets = []
        for b in options.bets.split(','):
            number = r'(\d+(?:\.\d*)?|\.\d+)'
            m = re.match('^%s$'%number,b)
            if m:
                bets.append(float(m.group(1)))
            else:
                m = re.match('^%s:%s:%s$'%(number,number,number),b)
                if m:
                    bets.extend(arange(float(m.group(1)),float(m.group(2)),float(m.group(3))))
                else:
                    print>>sys.stderr,"invalid bet sequence '%s'"%b
                    sys.exit(0)
        bets = array(bets)
    else:
        bets = bet/options.variable*arange(options.variable+1)
    print 'blinds = %g %g'%(sb,bb)
    print 'bets = %s'%bets
    cache = 'variable-sb%g-bb%g-bets%d-%d.npz'%(sb,bb,len(bets),abs(ahash(sb,bb,bets)))
    print 'cache file = %s'%cache
    print 'command = %s'%' '.join(sys.argv)
    print
    if os.path.exists(cache):
        data = load(cache)
        assert allclose(sb,data['sb'])
        assert allclose(bb,data['bb'])
        assert allclose(bets,data['bets'])
        equity,alice,bob = data['equity'],data['alice'],data['bob']
    else:
        equity,alice,bob = variable_bet_nash_equilibrium(sb,bb,bets)
        savez('variable-sb%g-bb%g-bets%d-%d.npz'%(sb,bb,len(bets),abs(ahash(sb,bb,bets))),sb=sb,bb=bb,bets=bets,equity=equity,alice=alice,bob=bob)
    _,alice_equity,bob_equity = check_variable_bet_nash_equilibrium(sb,bb,bets,equity,alice,bob)
    def times(p,s,tol=1e-3):
        if p>=1-tol:
            return s+' '
        if p<=tol:
            return ''
        return '%g*%s '%(p,s)
    print "\nAlice's strategy:"
    order = argsort(alice_equity.max(axis=0))[::-1]
    for j in order:
        print '  %s : equity %g, strategy %s%s'%(hands[j],alice_equity[:,j].max(),times(1-alice[:,j].sum(),'f'),''.join(times(alice[i,j],'b%g'%b) for i,b in enumerate(bets)))
    print "\nBob's strategy:"
    for i,b in enumerate(bets):
        print '  b%g : %s'%(b,''.join(times(bob[i,k],hands[k]) for k in order))
    if options.type:
        import pylab
        if options.type=='p':
            pb = sum(approx_hand_prob.reshape(1,-1)*alice,axis=1)
            pb = hstack([1-pb.sum(),pb])
            pylab.bar(hstack([-2,bets])-.5,pb)
            pylab.xlim(-1,bets.max())
            pylab.xlabel('action: raise amount or -1 for fold')
            pylab.ylabel('probability')
            pylab.title("Alice's probability of each action")
            pylab.show()
        else:
            raise RuntimeError("unknown plot type '%s'"%options.type)
else:
    equity,alice,bob = poker_nash_equilibrium(bet)
    print "Alice's exact payoff = %s (%g)"%(equity,equity)