[WIP] Add AS algorithm.

[shizejin]

An implementation of AS algorithm proposed by Abreu and Sannikov 2014. This algorithm is for computing the set V* of payoff pairs of all pure-strategy subgame-perfect equilibria with public randomization. Comparing to the JYC algorithm (which has been nicely implemented by @cc7768), this algorithm has significant efficiency gains, which arise from a better understanding of the manner in which extreme points of the equilibrium payoffs set are generated (in other word, this algorithm drops a lot of points that will not be included in V*, which will be kept by JYC algorithm during the iteration process).

Here demonstrates the usage of this algorithm with 3 examples. The first two are examples from the paper, and the third is a simple Prisoner Dilemma game that used in the notebook of JYC algorithm. I did comparison between JYC and AS algorithm using the third example.

I use Polyhedra.jl to do computation of polyhedra. Because this algorithm computes polyhedra intensively, the choice of PolyhedraLibrary will influence the speed very much. As the result shows, in the simple Prisoner Dilemma game, the AS algorithm is approximately twice faster than JYC algorithm if we use SimplePolyhedraLibrary, and 10 times faster if we use LRSLibrary.

However, the improvement in speed is way lower comparing to Java implementation in the paper. There are quite a few room for improving the efficiency. I am still working on this.

Detailed docstring and tests will be added in recent.

I really appreciate for any comments and suggestions.
@oyamad @QBatista @cc7768 @sglyon

[coveralls]

Coverage decreased (-12.4%) to 81.701% when pulling 9cad7ad on shizejin:add_AS_algorithm into 7ffab27 on QuantEcon:master.

[codecov-io]

Codecov Report

Merging #65 into master will increase coverage by 0.52%.
The diff coverage is 98.03%.

@@            Coverage Diff             @@
##           master      #65      +/-   ##
==========================================
+ Coverage   94.06%   94.58%   +0.52%     
==========================================
  Files           5        5              
  Lines         337      388      +51     
==========================================
+ Hits          317      367      +50     
- Misses         20       21       +1

Impacted Files	Coverage Δ
src/repeated_game.jl	93.33% <98.03%> (+2.1%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7ffab27...4fbf004. Read the comment docs.

[coveralls]

Coverage decreased (-12.4%) to 81.701% when pulling 9c71820 on shizejin:add_AS_algorithm into 7ffab27 on QuantEcon:master.

[cc7768]

Wow! This is great @shizejin -- When I wrote the original outerapproximation code, I had meant to do this algorithm as well but never got around to it. I am very excited that you have implemented it.

I've downloaded your notebook and can reproduce all of your graphs.

This is a first glance through the code.

Overall, I think things look very good. I'd love to see if we can make this code as fast as their implementation (or faster!).

Thanks @shizejin .

[cc7768]

If you don't know how many vertices you will need beforehand, I recommend using sizehint!(v_new, XYZ) where XYZ is maximum number of vertices you might have. If I remember correctly, they are able to pin down maximum number of vertices possible for finite action games (It might have been 4 * nactions?).

[shizejin]

Yes, the maximum number of vertices is 4 * num_actionprofiles. I will use sizehint!(v_new, 8 * num_actionprofiles) then. Thanks!

[cc7768]

These lines seem to be creating some sort of type instability -- I couldn't quite work out where it was coming from, but given that you call these functions again later, you might worry that the performance losses are coming from this instability (I'm not totally sure though)

julia> @code_warntype AS(rpd, 1000, plib, 1e-5, [0.0, 0.0])
Variables:
  #self#::Games.#AS
  rpd::Games.RepeatedGame{2,Float64}
  max_iter::Int64
  plib::Polyhedra.SimplePolyhedraLibrary{Float64}
  tol::Float64
  u@_6::Array{Float64,1}
  i@_7::Any
  #temp#@_8::Any
  a1::Int64
  payoff1::Float64
  payoff2::Float64
  IC1::Float64
  IC2::Float64
  p_IC::Polyhedra.SimplePolyhedron{_,Float64} where _
  p_inter::Any
  #temp#@_16::Int64
  a2::Int64
  #temp#@_18::Int64
  iter@_19::Int64
  v_new::Any
  u_::Array{Float64,2}
  #temp#@_22::Int64
  v_old::Array{Float64,2}
  H::Polyhedra.SimpleHRepresentation{_,Float64} where _
  V::Polyhedra.SimpleVRepresentation{_,_} where _ where _
  p::Polyhedra.SimplePolyhedron{_,Float64} where _
...

[shizejin]

This part is also what I am really concerned about. What I am trying to do here is to use a Matrix v_new to create VRepresentation without any redundant points (because the computation of intersection of two polyhedra is really slow if there are many redundant points).

I am not sure why Polyhedra.jl do not implement removevredundancy!(p::VRep), and I have to call the internal function removevredundancy(vrep::VRep, hrep::HRep; strongly=true) to remove the redundant points, which I don't know if it's efficient or not. I think we need to look into Polyhedra.jl deeply to see if I can improve this (very important) part.

[cc7768]

Could you precompute these?

Something like

p_IC = [polyhedron(SimpleHRepresentation(-eye(2), [-IC1[i], -IC2[i]]), plib) for i in naction]

where IC1 and IC2 were appropriately defined vectors? The IC constraints don't change, right? This might help with type stability since they'll be known before the for loop? I'm uncertain about this though

[shizejin]

IC constraints will change when u is updated. But yeah, I think there is no need to update p_IC in every iteration. I will move this into here so that p_IC is updated only if u changes. This will save a lot of time! Thanks.

[cc7768]

Given that u is always one dimensional, I'd put a comma here to make u_ one dimensional as well.

Below when you write u = u_ you are changing u between a 1 dimensional array and 2 dimensional array (this might also be a source of instability). See below for example

julia> u = [rand(), rand()]
2-element Array{Float64,1}:
 0.0212796
 0.144139 

julia> u_ = [rand() rand()]
1×2 Array{Float64,2}:
 0.8921  0.408161

julia> if all(u_ .> u)
           u = u_
       end
1×2 Array{Float64,2}:
 0.8921  0.408161
julia> u
1×2 Array{Float64,2}:
 0.8921  0.408161

[shizejin]

You are right. At the beginning I set both u and u_ to be two dimensional, but then I change u to one dimensional, while forget to change u_ as well. Thanks for finding this out!

[shizejin]

Thanks for your review @cc7768! It is very helpful.

There are two things in my mind that might improve the efficiency:

1. The way to get rid of redundant points. Now I am first calculating all the points, and then use removevredundancy to remove the redundant points. Maybe another way is to construct a new VRepresentation at the beginning of each iteration, and update it gradually. (If the found point is in the corresponding HRepresentation, we do not change the VRepresentation. Otherwise, we add this point to VRepresentation. One possible drawback is that we have to transform VRepresentation to HRepresentation for many times.) I am not sure which way would be faster.

2. The way to find intersection points for each pair of IC boundaries. Now I am constructing a polyhedron p_IC and then use method intersect. After get the vertices of intersection polyhedron p_inter, I check each vertex to see if v[1] == IC1 or v[2] == IC2. Another possible way is to set a optimization problem. Given the current polyhedron p, what is the extreme value of v[2] if we set v[1] = IC1. In this case, we don't need to construct p_IC and don't need to intersect polyhedra neither. I am not sure if Polyhedra.jl provides method to do this though.

[coveralls]

Coverage increased (+0.5%) to 94.588% when pulling 7d4b2ee on shizejin:add_AS_algorithm into 7ffab27 on QuantEcon:master.

[coveralls]

Coverage increased (+0.5%) to 94.588% when pulling 4fbf004 on shizejin:add_AS_algorithm into 7ffab27 on QuantEcon:master.

[sglyon]

Thank you @shizejin for writing this up.

Also thanks to @cc7768 for the initial review. @cc7768, would you be willing to see this PR through to the end and eventually merge when it is ready to go?

[oyamad]

@shizejin What is the status of this PR?

[shizejin]

@oyamad Currently the performance is still not as desirable. But it seems that there has been many updates in Polyhedra.jl. I will check in a few days to see if these updates solve the performance issue in this PR.