# naive score function; discrete points {0, 0.5, 1}; unrealistic 
function get_simple_score(player, opponent, order, round, elo)
    x = elo[player] - elo[opponent] + 0.75 * (-1)^(order + 1)
    if x > 0 
        return 1
    elseif x < 0
        return 0
    else 
        return 0.5
    end  
end;

# expected score of a game (continuous); each game identical; first-move advantage;
function get_chess_score(player, opponent, order, round, elo)
    x = elo[player] - elo[opponent] + 35*(-1)^(order + 1)
    return 1/(1 + 10^(-x/400)) 
end;

# expected score of a game; each game being distinct; first-move advantage;
function get_chess_score2(player, opponent, order, round, elo)
    x = elo[player] - elo[opponent] + 35*(-1)^(order % 2 + 1)      
    return 1/(1 + 10^(-x/400)) * (1.01 - 0.01round/N_PLAYERS)
end;

#  expected score of a game; each game being distinct; selection from a pool of players 
function get_chess_score3(player, opponent, order, round, elo)
    x = elo[player] - elo[opponent] + 35*(-1)^(order % 2 + 1)      
    return 1/(1 + 10^(-x/400)) * (1.01 - 0.01round/N_POOL)
end;

# total expected score of match; each game ientical
function get_match_score(team_player, team_opponent, order, elo)
    sum(get_chess_score(team_player[i], team_opponent[i], 
            i + order + 1, elo) for i=1:length(team_player))
end;
 
# total expected score of match; each game distinct
function get_match_score2(team_player, team_opponent, order, elo)
    sum(get_chess_score2(team_player[i], team_opponent[i], 
            i + order + 1, i, elo) for i=1:length(team_player))
end;

# one is enough!! No of games - 1's score = 2's score
#Pkg.add("Images") # add images package if not already installed
using JuMP, Clp, Cbc, NamedArrays
using PyPlot
import Combinatorics

raw = readcsv("./data/Chess.csv")
(m,n) = size(raw) # dimensions of data

players = raw[2:end, 2]      # vector of players
elos = raw[2:end, 3]         # vector of Elo ratings
federations = raw[2:end, 6]  # vector of 
worth = raw[2:end, 9]        # player worth in thousand dolaars
num_players = length(players)# number of total players (200)

player_elo = Dict(zip(players, elos))           # Dictionary player to Elo rating
player_country = Dict(zip(players, federations))# Dictionary player to country
player_worth = Dict(zip(players, worth))        # Dictionary player to worth

player_pool_USA = []
player_pool_Russia = []

# building USA and Russia pools
for i in 2: m
    if raw[i, 6] == "United States"
        push!(player_pool_USA, raw[i, 2])
    elseif raw[i, 6] == "Russia"
        push!(player_pool_Russia, raw[i, 2])
    end
end

N_PLAYERS = 4 # number of players in both teams
team_USA = player_pool_USA[1:N_PLAYERS]       # Best 4 USA players form the home team
team_Russia = player_pool_Russia[1:N_PLAYERS];# Best 4 Russia players form the visitor team

# sticking to descending elo order for fixed strategy
plan_Russia = [1 2 3 4];

pool = players
# N_POOL = length(pool) # top 200 players
N_POOL = 10    # top 10 players
N_CHOSEN = 4;  # teams of 4 to be formed

function print_match(aopt::Array, team_name, players, int_mode=true)
    m = size(aopt)[1]
    n = size(aopt)[2]
    if int_mode
        solution = NamedArray([Int(aopt[i, k]) 
                for i=1:m, k=1:n], (players, [1:n;]),(team_name, "Table"))
    else
        solution = NamedArray([aopt[i, k]
                for i=1:m, k=1:n], (players, [1:n;]),(team_name, "Table"))
    end
    println(solution)
    println()
end

function print_fixed_strategy(xopt, N_PLAYERS, 
        team_home, team_visitor, plan_visitor, elo)   
    println("Strategy: ")
    for k in 1:N_PLAYERS
        for i in 1:N_PLAYERS
            if xopt[i, k] != 0
                @printf "%20s (%d) vs " team_home[i] elo[team_home[i]] 
                @printf "(%d) %s\n" elo[team_visitor[plan_visitor[k]]] team_visitor[plan_visitor[k]]
            end
        end
    end
end
;

function print_pure_strategy(N_PLAYERS, team_home, team_visitor, elo, xopt, yopt)
    for k in 1:N_PLAYERS # round
        print("Round: ", k)
        for i in 1:N_PLAYERS
            if xopt[i, k] != 0
                home_player = string(team_home[i], " (",elo[team_home[i]], ")")                
                @printf "%30s %5.5s " home_player "vs."
            end
        end
        for j in 1:N_PLAYERS
            if yopt[j, k] != 0
                visitor_player = string( "(",elo[team_visitor[j]], ") ", team_visitor[j])
                @printf " %s\n" visitor_player
            end
        end        
    end
end;

# score_matrix function from 2 lists
function get_score_matrix(list1::Array, list2::Array, players_dict, 
        get_match)
    # assume that team 1 plays first
    dim1 = length(list1)
    dim2 = length(list2)
    scores = Matrix(dim1, dim2)
    for i in 1:dim1
        for j in 1:dim2
            scores[i, j] = get_match(list1[i], list2[j], 1, players_dict)
        end
    end
    return scores
end;

# ordered combination: return a list of nplayers obtained from a pool
function get_combinations(list_players::Array, nchosen::Int)
    return collect(combinations(list_players, nchosen))
end;

function get_permutations(list_players::Array)
    return collect(permutations(list_players, N_PLAYERS))
end;

list_home = get_permutations(team_USA[1:N_PLAYERS])
list_visitor = get_permutations(team_Russia[1:N_PLAYERS])

scores_matrix = get_score_matrix(list_home, list_visitor, player_elo, get_match_score2);

function get_fixed_strategy(N_PLAYERS, team_home, 
        team_visitor, plan_visitor, get_score, elo)
    simpleModel = Model(solver=CbcSolver())

    @variable(simpleModel, x[1:N_PLAYERS, 1:N_PLAYERS], Bin) # matrix of main assignment variables

    @constraint(simpleModel, supply[k in 1:N_PLAYERS], 
                    sum(x[k, j] for j=1:N_PLAYERS) == 1)    # supply constraint
    @constraint(simpleModel, demand[k in 1:N_PLAYERS], 
                    sum(x[i, k] for i=1:N_PLAYERS) == 1)    # demand constraint

    @objective(simpleModel, Max, sum(x[i, k] * get_score(
                team_home[i], team_visitor[plan_visitor[k]], k, k, elo) 
                for i=1:N_PLAYERS, k=1:N_PLAYERS))          # maximization objective function

    solve(simpleModel)
    println("Score: ", getobjectivevalue(simpleModel))
    xopt = getvalue(x)
    return xopt
end;

xopt = get_fixed_strategy(N_PLAYERS, 
    team_USA, team_Russia, plan_Russia, get_simple_score, player_elo);
print_fixed_strategy(xopt, N_PLAYERS, 
    team_USA, team_Russia, plan_Russia, player_elo)

xopt = get_fixed_strategy(N_PLAYERS, 
    team_USA, team_Russia, plan_Russia, get_chess_score, player_elo);
print_fixed_strategy(xopt, N_PLAYERS, 
    team_USA, team_Russia, plan_Russia, player_elo)

x_opt = get_fixed_strategy(N_PLAYERS,
    team_USA, team_Russia,  plan_Russia, get_chess_score2, player_elo);
print_fixed_strategy(xopt, N_PLAYERS, 
    team_USA, team_Russia, plan_Russia, player_elo)

function get_team_pure_strategy_maximin(N_PLAYERS, team_home, team_visitor, get_score, elo)
    maximinModel = Model(solver=CbcSolver())

    @variable(maximinModel, x[1:N_PLAYERS, 1:N_PLAYERS], Bin)  # matrix of main assignment variables
    @variable(maximinModel, p[1:N_PLAYERS])                    # vector of inner model dual variables
    @variable(maximinModel, q[1:N_PLAYERS])                    # vector of inner model dual variables
    @variable(maximinModel, r[1:N_PLAYERS, 1:N_PLAYERS] <= 0)  # vector of inner model dual variables

    @constraint(maximinModel, supply[k in 1:N_PLAYERS], 
                    sum(x[k, j] for j=1:N_PLAYERS) == 1)       # supply constraint
    @constraint(maximinModel, demand[k in 1:N_PLAYERS], 
                    sum(x[i, k] for i=1:N_PLAYERS) == 1)       # demand constraint 

    @constraint(maximinModel, mincons[j in 1:N_PLAYERS, k in 1:N_PLAYERS], 
                p[j] + q[k] + r[j, k] <= 
                sum(get_score(team_home[i], team_visitor[j], k, k, elo) * x[i, k] 
                for i=1:N_PLAYERS))                            # inner dual constraint

    @objective(maximinModel, Max, sum(p) + sum(q) + sum(r))    # maximization objective function    

    solve(maximinModel)
                                        
    println("Maximin Score: ", getobjectivevalue(maximinModel))
                                        
    return getvalue(x)                                    
end;

function get_expected_opponent_pure_strategy_maximin(N_PLAYERS, team_home, team_visitor, get_score, elo, xopt)
    ym = Model(solver=CbcSolver())

    # matrix of assignment variables
    @variable(ym, y[1:N_PLAYERS, 1:N_PLAYERS], Bin)

    # supply and demand constraints
    @constraint(ym, supply[k in 1:N_PLAYERS], sum(y[k, j] for j=1:N_PLAYERS) == 1)
    @constraint(ym, demand[k in 1:N_PLAYERS], sum(y[i, k] for i=1:N_PLAYERS) == 1)

    # objective expression
    @expression(ym, score, sum(get_score(team_home[i], 
                        team_visitor[j], k, k, elo) * xopt[i, k] * y[j, k] 
                        for i=1:N_PLAYERS, j in 1:N_PLAYERS, k in 1:N_PLAYERS))
    # minimization of objective expression
    @objective(ym, Min, score)

    solve(ym)

    return getvalue(y)
end;

function get_opponent_pure_strategy_minimax(N_PLAYERS, team_home, 
        team_visitor, get_score, elo)
    
    minimaxModel = Model(solver=CbcSolver())

    @variable(minimaxModel, y[1:N_PLAYERS, 1:N_PLAYERS], Bin) # matrix of main assignment variables
    @variable(minimaxModel, p[1:N_PLAYERS])                   # vector of inner model dual variables
    @variable(minimaxModel, q[1:N_PLAYERS])                   # vector of inner model dual variables
    @variable(minimaxModel, r[1:N_PLAYERS, 1:N_PLAYERS] >= 0) # vector of inner model dual variables

    @constraint(minimaxModel, supply[k in 1:N_PLAYERS], 
                    sum(y[k, j] for j=1:N_PLAYERS) == 1)      # supply constraint
    @constraint(minimaxModel, demand[k in 1:N_PLAYERS], 
                    sum(y[i, k] for i=1:N_PLAYERS) == 1)      # demand constraint

    @constraint(minimaxModel, mincons[j in 1:N_PLAYERS, k in 1:N_PLAYERS], 
                p[j] + q[k] + r[j, k] >= 
                sum(get_score(team_home[i], team_visitor[j], k, k, elo) * y[i, k] 
                for i=1:N_PLAYERS))                           # inner dual constraint

    @objective(minimaxModel, Min, sum(p) + sum(q) + sum(r))   # minimization objective function

    solve(minimaxModel)

    println("Minimax Score: ", getobjectivevalue(minimaxModel))
    
    return getvalue(y)
end;
                                    
function get_expected_team_pure_strategy_minimax(N_PLAYERS, team_home, team_visitor, get_score, elo, yopt)
    xm = Model(solver=CbcSolver())
    # assignment variable
    @variable(xm, x[1:N_PLAYERS, 1:N_PLAYERS], Bin)
    # supply and demand constraints
    @constraint(xm, supply[k in 1:N_PLAYERS], sum(x[k, j] for j=1:N_PLAYERS) == 1)
    @constraint(xm, demand[k in 1:N_PLAYERS], sum(x[i, k] for i=1:N_PLAYERS) == 1)
    # objective expression
    @expression(xm, score, sum(get_score(team_home[i], 
                        team_visitor[j], k, k, elo) * yopt[i, k] * x[j, k] 
                        for i=1:N_PLAYERS, j in 1:N_PLAYERS, k in 1:N_PLAYERS))
    # maximization of objective expression
    @objective(xm, Max, score)

    solve(xm)

    return getvalue(x)
end;

xopt = get_team_pure_strategy_maximin(N_PLAYERS, 
    team_USA, team_Russia, get_chess_score, player_elo);
println("Home Strategy")
print_match(xopt, "USA", team_USA)

yopt = get_expected_opponent_pure_strategy_maximin(N_PLAYERS, 
    team_USA, team_Russia, get_chess_score, player_elo, xopt)
println("Expected Visitor Strategy")
print_match(yopt, "Russia", team_Russia)

println("Maximin Strategy: ")
print_pure_strategy(N_PLAYERS, team_USA, team_Russia, player_elo, xopt, yopt)  

yopt = get_opponent_pure_strategy_minimax(N_PLAYERS, 
    team_USA, team_Russia, get_chess_score, player_elo);

println("Visitor Strategy")
print_match(yopt, "Russia", team_Russia)

xopt =  get_expected_team_pure_strategy_minimax(N_PLAYERS, 
    team_USA, team_Russia, get_chess_score, player_elo, yopt);
println("Expected Home Strategy")
print_match(xopt, "USA", team_USA) 

print_pure_strategy(N_PLAYERS, team_USA, team_Russia, player_elo, xopt, yopt)

function plot_heatmap(scores_matrix::Array, labelx, labely)
    imshow(scores_matrix, cmap="jet", aspect="auto", interpolation="nearest")
    colorbar()
    plt[:ylabel](labely)
    plt[:xlabel](labelx)
end;
plot_heatmap(scores_matrix, "Russia", "USA");

function get_MSNE(N_PLAYERS, team_home, team_visitor, get_score, elo)
    nashModel = Model(solver=ClpSolver())

    @variable(nashModel, 1 >= x[1:N_PLAYERS, 1:N_PLAYERS] >= 0)# matrix of main assignment variables
    @variable(nashModel, p[1:N_PLAYERS])                       # vector of inner model dual variables 
    @variable(nashModel, q[1:N_PLAYERS])                       # vector of inner model dual variables
    @variable(nashModel, r[1:N_PLAYERS, 1:N_PLAYERS] <= 0)     # vector of inner model dual variables

    @constraint(nashModel, supply[k in 1:N_PLAYERS], 
                    sum(x[k, j] for j=1:N_PLAYERS) == 1)       # supply constraint
    @constraint(nashModel, demand[k in 1:N_PLAYERS], 
                    sum(x[i, k] for i=1:N_PLAYERS) == 1)       # supply constraint 

    @constraint(nashModel, mincons[j in 1:N_PLAYERS, k in 1:N_PLAYERS], 
                p[j] + q[k] + r[j, k] <= 
                sum(get_score(team_home[i], team_visitor[j], k, k, elo) 
                    * x[i, k] for i=1:N_PLAYERS))              # inner model dual constraint

    @objective(nashModel, Max, sum(p) + sum(q) + sum(r))       # maximization objective function     

    solve(nashModel)

    println("Optimal Score: ", getobjectivevalue(nashModel))
                                        
    return getvalue(x)
end;

function get_MSNE_opponent(N_PLAYERS, team_home, team_visitor, get_score, elo)
    nashModel = Model(solver=ClpSolver())

    @variable(nashModel, 0 <= y[1:N_PLAYERS, 1:N_PLAYERS] <= 1)# matrix of main assignment variables
    @variable(nashModel, p[1:N_PLAYERS])                       # vector of inner model dual variables
    @variable(nashModel, q[1:N_PLAYERS])                       # vector of inner model dual variables
    @variable(nashModel, r[1:N_PLAYERS, 1:N_PLAYERS] >= 0)     # vector of inner model dual variables

    @constraint(nashModel, supply[k in 1:N_PLAYERS], 
                    sum(y[k, j] for j=1:N_PLAYERS) == 1)       # supply constraint
    @constraint(nashModel, demand[k in 1:N_PLAYERS], 
                    sum(y[i, k] for i=1:N_PLAYERS) == 1)       # demand constraint

    @constraint(nashModel, maxcons[i in 1:N_PLAYERS, k in 1:N_PLAYERS],
                p[i] + q[k] + r[i, k] >= 
                sum(get_score(team_home[i], team_visitor[j], k%2, k, elo) 
                    * y[j, k] for j=1:N_PLAYERS))              # inner model dual constraint

    @objective(nashModel, Min, sum(p) + sum(q) + sum(r))       # minimization objective function      

    solve(nashModel)

    println("Optimal Score: ", getobjectivevalue(nashModel))
                                        
    return getvalue(y)
end;

xopt = get_MSNE(N_PLAYERS, 
    team_USA, team_Russia, get_chess_score2, player_elo);
println("Home Mixed Strategy")
print_match(xopt, "USA", team_USA, false)

yopt = get_MSNE_opponent(N_PLAYERS, 
    team_USA, team_Russia, get_chess_score2, player_elo);
println("Visitor Mixed Strategy")
print_match(yopt, "Russia", team_Russia, false)

### Select a team compsition from a pool
function select_from_pool_mixed(N_POOL, N_CHOSEN, get_score, elo)
    inf = -1000
    sup = 1000 #

    xm = Model(solver=CbcSolver())

    # matrix of main probabilistic assignment variables
    @variable(xm, 0 <= x[1:N_POOL, 1:N_CHOSEN] <= 1)
    # matrix of main assignment variables
    @variable(xm, z[1:N_POOL, 1:N_CHOSEN], Bin)
    # vector of inner model dual variables
    @variable(xm, p[1:N_POOL] <= 0) 
    # vector of inner model dual variables
    @variable(xm, q[1:N_CHOSEN])
    # matrix of inner model dual variables
    @variable(xm, r[1:N_POOL, 1:N_CHOSEN] <= 0)
    # vector of epigraph variables
    @variable(xm, t[1:N_POOL] <= 0)

    # supply constraint of z
    @constraint(xm, supply[i in 1:N_POOL], 
                    sum(z[i, k] for k=1:N_CHOSEN) <= 1)
    # demand constraint of z
    @constraint(xm, demand[k in 1:N_CHOSEN], 
                    sum(z[i, k] for i=1:N_POOL) == 1)
    # supply constraint of x                   
    @constraint(xm, supplyProb[i in 1:N_POOL], 
                                sum(x[i, k] for k=1:N_CHOSEN) <= 1)
    # demand constraint of x
    @constraint(xm, demandProb[k in 1:N_CHOSEN], 
                                sum(x[i, k] for i=1:N_POOL) == 1)


    # z = 0 => x = 0
    @constraint(xm, node[i in 1:N_POOL], 
            sum(x[i, k] for k=1:N_CHOSEN) <= sum(z[i, k] for k=1:N_CHOSEN))
    # epigraph constraint 1                                                               
    @constraint(xm, icons[j in 1:N_POOL], t[j] >= 
                                (1 - sum(z[j, k] for k=1:N_CHOSEN))*inf)
    # epigraph constraint 2   
    @constraint(xm, ucons[j in 1:N_POOL], t[j] <= 
                                p[j] + sup*sum(z[j, k] for k=1:N_CHOSEN))
    # dual constraint
    @constraint(xm, mincons[j in 1:N_POOL, k in 1:N_CHOSEN], 
        p[j] + q[k] + r[j, k]  <= sum(get_score(i, j, k, k, elo) 
                                                    * x[i, k] for i=1:N_POOL))
    # objective
    @objective(xm, Max, sum(t) + sum(q) + sum(r))        

    solve(xm)
    println("Optimal Score: ", getobjectivevalue(xm))
    return getvalue(x)
end;

### Select a team compsition from a pool with a budget trade-off
function select_with_budget_mixed(lambda, N_POOL, N_CHOSEN, 
        get_score, elo, worth)
    inf = -1000
    sup = 1000 

    xm = Model(solver=CbcSolver())

    # matrix of main probabilistic assignment variables
    @variable(xm, 0 <= x[1:N_POOL, 1:N_CHOSEN] <= 1)
    # matrix of main assignment variables
    @variable(xm, z[1:N_POOL, 1:N_CHOSEN], Bin)
    # vector of inner model dual variables
    @variable(xm, p[1:N_POOL] <= 0) 
    # vector of inner model dual variables
    @variable(xm, q[1:N_CHOSEN])
    # matrix of inner model dual variables
    @variable(xm, r[1:N_POOL, 1:N_CHOSEN] <= 0)
    # vector of epigraph variables
    @variable(xm, t[1:N_POOL] <= 0)

    # supply constraint of z
    @constraint(xm, supply[i in 1:N_POOL], 
                    sum(z[i, k] for k=1:N_CHOSEN) <= 1)
    # demand constraint of z
    @constraint(xm, demand[k in 1:N_CHOSEN], 
                    sum(z[i, k] for i=1:N_POOL) == 1)
    # supply constraint of x                   
    @constraint(xm, supplyProb[i in 1:N_POOL], 
                                sum(x[i, k] for k=1:N_CHOSEN) <= 1)
    # demand constraint of x
    @constraint(xm, demandProb[k in 1:N_CHOSEN], 
                                sum(x[i, k] for i=1:N_POOL) == 1)

    # z = 0 => x = 0
    @constraint(xm, node[i in 1:N_POOL], 
            sum(x[i, k] for k=1:N_CHOSEN) <= sum(z[i, k] for k=1:N_CHOSEN))
                                                                    
    # epigraph constraint 1                                                               
    @constraint(xm, icons[j in 1:N_POOL], t[j] >= 
                                (1 - sum(z[j, k] for k=1:N_CHOSEN))*inf)
    # epigraph constraint 2   
    @constraint(xm, ucons[j in 1:N_POOL], t[j] <= 
                                p[j] + sup*sum(z[j, k] for k=1:N_CHOSEN))
    # dual constraint
    @constraint(xm, mincons[j in 1:N_POOL, k in 1:N_CHOSEN], 
        p[j] + q[k] + r[j, k]  <= sum(get_score(i, j, k, k, elo) 
                                                    * x[i, k] for i=1:N_POOL))
    # Player cost expression
    @expression(xm, playercost,  lambda * 
                sum(sum(z[i, k] for k=1:N_CHOSEN) 
                        * worth[i] for i=1:N_POOL))
    # Objective function
    @objective(xm, Max, sum(t) + sum(q) + sum(r) - playercost)        

    solve(xm)

#     println("Optimal Score: ", getobjectivevalue(xm))
    return [getvalue(x), getobjectivevalue(xm), getvalue(playercost)] 
end;

xopt = select_from_pool_mixed(N_POOL, N_CHOSEN, get_chess_score3, elos);
print_match(xopt, "Name", players[1:N_POOL], false)

function explore_tradeoff(lambda)   
    result = select_with_budget_mixed(lambda, N_POOL, N_CHOSEN, 
        get_chess_score3, elos, worth * 0.001)
    xopt = result[1]
    cost = result[3]
    score = result[2] + cost
    println("Expected score: ", score)
    println("Cost: ", cost)
    print_match(xopt, "Name", players[1:N_POOL], false)
end;

explore_tradeoff(0.1) 

explore_tradeoff(0.2) 

explore_tradeoff(1)

# compute optimal tradeoff curve (this may take a few seconds
N = 30
lambda_values = linspace(0.1,1,N)
# lambda_values = [0.001, 0.05, 0.01, 0.1, 1, 10]
N = length(lambda_values)
cost = zeros(N)
score = zeros(N)
opt = zeros(N)
for (i, lambda) in enumerate(lambda_values)
   result = select_with_budget_mixed(lambda, N_POOL, N_CHOSEN, 
        get_chess_score3, elos, worth * 0.001)
    cost[i] = result[3]
    score[i] = result[2] + result[3]
    opt[i] = result[2]
end

using PyPlot
plot(cost, score,"b-")
xlabel("cost (k)")
ylabel("expected score")
title("Pareto Curve")
tight_layout()

title("Expected Score vs. Lambda")
plot(lambda_values, score,"g-")
xlabel("lambda")
ylabel("expected score");

function get_maximin_pure_by_enumeration(score_matrix::Array, maxPlayer=true)
    n = size(scores_matrix)[1]
    t = zeros(n)
    if maxPlayer             
        for k in 1:n
            t[k] = findmin(score_matrix[k, 1:n])[1]
        end
        return findmax(t)
    else
        for k in 1:n
            t[k] = findmax(score_matrix[1:n, k])[1]
        end
        return findmin(t)
    end
end
;
println("Maximin Score: ",get_maximin_pure_by_enumeration(scores_matrix, true)) # true for primal
println("Minimax Score: ",get_maximin_pure_by_enumeration(scores_matrix, false)) # false for primal

function get_MSNE_by_enumeration(score_matrix::Array)
    n = size(scores_matrix)[1]
    # for primal - P person
    m = Model(solver=ClpSolver())

    @variable(m, 1>= p[1:n] >=0)
    @variable(m, t)

    @constraint(m, sum(p) == 1)
    @constraint(m, min_const[k in 1:n], 
        sum(p[i] * score_matrix[i, k] for i=1:n) >= t)
    # Max min (g1, g2, ..., gn)
    @objective(m, Max, t)
    solve(m)
                
    popt = getvalue(p)
#     println(popt)
    println(getobjectivevalue(m))
    return popt
end;
function get_MSNE_dual_by_enumeration(score_matrix::Array)
    n = size(scores_matrix)[2]
    # for primal - P person
    m = Model(solver=ClpSolver())

    @variable(m, 1>= q[1:n] >=0)
    @variable(m, t)

    @constraint(m, sum(q) == 1)
    @constraint(m, max_const[k in 1:n], 
        sum(q[i] * score_matrix[k, i] for i=1:n) <= t)
    # Max min (g1, g2, ..., gn)
    @objective(m, Min, t)
    solve(m)
                
    qopt = getvalue(q)
#     println(qopt)
    println(getobjectivevalue(m))
    return qopt
end;  
        
print("Primal optimal score: ") 
popt = get_MSNE_by_enumeration(scores_matrix);
print("Dual optimal score: ") 
popt =get_MSNE_dual_by_enumeration(scores_matrix);;