#Pkg.clone("https://github.com/squipbar/Polygons.jl")
#Pkg.clone("https://github.com/cc7768/CHull2d.jl")

using Polygons

Z = [ 1 -1; 1 1; -1 1; -1 -1 ]
G = [ 1 0; 0 1; -1 0; 0 -1 ]
m = [ 1, 1, 1, 1 ]
a = Polygon( Z, G, m ) ;
b = Polygon( pts=Z ) ;
c = Polygon( dirs=G, dists=m ) ;

badZ = Z[ [1, 3, 4, 2], : ]
d = Polygon(pts=badZ)

a

polyPlot(a)

f = a + [ 1, 1 ] ;

g = 1.5 * a ;

polyPlot( [a, f, g] )

h = Polygon( pts=[ 0 1 ; 1 0 ; -1 0 ; 0 -1 ] )
exactdirs = [ a.dirs ; h.dirs ]
# The union of the face directions.  Will give an exact result
exact = setSum( a, h, exactdirs )
# The exact set sum
apxdirs = [ .1 1; .2 .8 ; -1 1 ; -.7 -.3 ; .4 -.6; .6 -.4  ]
outer = setSum( a, h, apxdirs ) ;
inner = setSum( a, h, apxdirs, false ) ;
polyPlot( [ a, h, exact, inner, outer ])

N = 40
theta = 2 * pi * rand( N )
manydirs = zeros( N, 2 )
for i in 1:N
    manydirs[i,:] = [ cos(theta[i]), sin(theta[i]) ]
end
manyouter = setSum( a, h, manydirs ) ;
manyinner = setSum( a, h, manydirs, false ) ;
polyPlot( [ exact, manyinner, manyouter ])

using Gadfly
using Colors
N = 40
srand(42)
manypts = randn( N, 2 )
chull( manypts )
ch = Polygon( pts = chull( manypts ) )
chplot = ch.pts[ [ 1:end; 1], : ]
plot( layer( x=manypts[:,1], y=manypts[:,2], Geom.point ), layer( x=chplot[:,1], y=chplot[:,2],  Geom.path,
          Theme(default_color=colorant"red") ) ) 

chr = crop( ch, 1, .5, true )
chdown = crop( ch, 2, .5, false )
polyPlot( [chr, chdown, ch])

polyPlot( [ exact, ch, union( [ exact, ch ] ) ] )

P = [ .8 .2; .2 .8 ]
U = [ .5 -.5 ; -.5 .5 ]
bet = 0.9
vlow = ( 1 - bet ) * ( ( eye(2) - bet * P ) \ U )

N=40
dirs = zeros( N, 2 )
dirs[:,1] = [ cos(2 * pi * (i-1)/N) for i in 1:N ]
dirs[:,2] = [ sin(2 * pi * (i-1)/N) for i in 1:N ]
# Search directions.  Need to be careful to make sure that this is a matrix type.
U1(s) = [ 2 -1; 3+s s ]
U2(s) = [ 2 3-s; -1 -s ]
# The period payoff functions

function T_unc( W::Array{ Polygon, 1 }, s::Float64, a::Vector, outer = true)
# s is the value for the state and a is a vector containing the action indices
    u = [ U1(s)[a[1],a[2]] , U2(s)[a[1],a[2]] ] 
        # Period payoff
    s_idx = ( s == .5 ) ? 1 : 2
        # The index of s.  Required for selecting the right transition probability row
    return ( ( 1 - bet ) * u ) + ( bet * ( weightedSum( W, vec( P[ s_idx, :] ), dirs, outer ) ) )
end ;

winit1 = Polygon( pts = [ 2.0 2.0; -1.0 3.5; 3.5 -1.0; -0.5 0.5 ; 0.5 -0.5 ] ) 
winit = [ winit1, winit1 ]
uncouter = T_unc( winit, .5, [1, 1] )
uncinner = T_unc( winit, .5, [1, 1], false )
polyPlot( [winit1, uncouter, uncinner ])

function T_a( W::Array{ Polygon, 1 }, s::Float64, a::Vector, vlow::Matrix, outer = true)
    # calcultes the incentive-compatible value set conditional on the action a
    unc = T_unc( W, s, a, outer )
        # The unconstrained set
    pddev = [ maximum( U1(s)[ :, a[2] ] ), maximum( U2(s)[ a[1], : ] ) ]
        # The deviating period payoffs
    s_idx = ( s == .5 ) ? 1 : 2
        # The index of s.  Required for selecting the right transition probability row
    dev = ( 1 - bet ) * pddev + bet * vlow[ s_idx, : ]'
        # The deviating payoff
    ic = crop( crop( unc, 1, dev[1] ), 2, dev[2] )
        # The incentive compatible set of payoffs\
    return ic
end ;

icouter = T_a( winit, .5, [1, 1], vlow )
icinner = T_a( winit, .5, [1, 1], vlow, false )
polyPlot( [icouter, icinner ])

A_idx = [ 1 1; 1 2 ; 2 1; 2 2 ]
    # The matrix of action combinations

function T_operator( W::Array{ Polygon, 1 }, states::Vector, vlow, outer = true)
    # Forms the T operator
    out = [ union( [ T_a( W, s, vec(A_idx[j,:]), vlow, outer)::Polygon for j in 1:size(A_idx)[1] ] )::Polygon for s in states ]
    return out
end

polyPlot( [ winit1; T_operator( winit, [ .5 , -.5 ], vlow ) ] )

niter = 40
N=100
dirs = zeros( N, 2 )
dirs[:,1] = [ cos(2 * pi * (i-1)/N) for i in 1:N ]
dirs[:,2] = [ sin(2 * pi * (i-1)/N) for i in 1:N ]

function eqm( W::Array{Polygon,1}, states::Vector, outer::Bool=true )
    U = [ states' ; -states' ]
    vlow = ( 1 - bet ) * ( ( eye(2) - bet * P ) \ U )
    for ( i in 1:niter )
        Wnew = T_operator( W, states, vlow, outer )
        W = Wnew
    end
    return W
end

@time wouter = eqm( winit, [ .5, -.5 ] )
@time winner = eqm( winit, [ .5, -.5 ], false )

polyPlot( [ wouter; winner ] )