text:
using NLPModelsJuMP
using JuMP

function Jason1()
  P1 = JuMP.Model()
  JuMP.@variable(P1, x[1:2], start = 1.0)
  JuMP.@constraint(P1, x[1]^2 + x[2]^2 - 1 <= 0)
  JuMP.@constraint(P1, x[1] >= 0)
  JuMP.@constraint(P1, x[2] >= 0)
  JuMP.@NLobjective(P1, Min, 0.5 * (x[1] - 2 / 3)^2)
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[1:2], start = 1.0)
  JuMP.@constraint(P2, x[2] + x[1] - 1 <= 0)
  JuMP.@constraint(P2, x[1] >= 0)
  JuMP.@constraint(P2, x[2] >= 0)
  JuMP.@NLobjective(P2, Min, 0.5 * (x[2] - 1 / 3)^2)
  Player2 = MathProgNLPModel(P2)
  sol = x -> (x[1] == 2 / 3 && x[2] == 1 / 3)

  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function Jason2()
  P1 = JuMP.Model()
  JuMP.@variable(P1, x[1:3], start = 1.0)
  JuMP.@constraint(P1, x[1] + x[2] + x[3] - 3 <= 0)
  JuMP.@constraint(P1, x[1] >= 0)
  JuMP.@constraint(P1, x[2] >= 0)
  JuMP.@constraint(P1, x[3] >= 0)
  JuMP.@NLobjective(P1, Min, 0.5 * (x[1] - 1)^2 + 0.5 * (x[2] - 1)^2)
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[1:3], start = 1.0)
  JuMP.@constraint(P2, -x[1] - x[2] - x[3] + 3 <= 0)
  JuMP.@constraint(P2, x[1]^2 + x[2]^2 + x[3]^2 - 9 <= 0)
  JuMP.@constraint(P2, x[1] >= 0)
  JuMP.@constraint(P2, x[2] >= 0)
  JuMP.@constraint(P2, x[3] >= 0)
  JuMP.@NLobjective(P2, Min, 0.5 * x[1] * x[2] * x[3]^2)
  Player2 = MathProgNLPModel(P2)
  sol = x -> (t = 1 - x[1]; (x[2] == 1 - t) && (x[3] == 1 + 2 * t))

  return GNEP(2, [2, 1], [Player1, Player2], sol = sol)
end

function Tangi1()
  ul(i) = [0.1; 0.1; 0.1][i]
  ux(i) = [Inf; Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + x[2] + x[3] <= 3)
  JuMP.@NLobjective(P1, Min, 0.5 * (x[1] - 1)^2 + 0.5 * (x[2] - 1)^2)
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@NLconstraint(P2, -x[1]^2 * x[3]^2 <= -0.5)
  JuMP.@NLconstraint(P2, -x[2]^2 * x[3]^2 <= -0.5)
  JuMP.@NLobjective(P2, Min, 0.5 * x[1] * x[2] * x[3]^2)
  Player2 = MathProgNLPModel(P2)

  sol = x -> (t = x[1]; (x[1] == 5 && x[2] == 9) || (x[2] == 15 - t && t <= 10 && t >= 9))
  return GNEP(2, [2, 1], [Player1, Player2], sol = sol)
end

function Jason2cvx()
  ul(i) = [log(0.1); log(0.1); log(0.1)][i]
  ux(i) = [Inf; Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, y[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  #JuMP.@NLconstraint(P1,exp(y[1]+log(1/3))+exp(y[2]+log(1/3))+exp(y[3]+log(1/3))<=1)
  #JuMP.@NLconstraint(P1,log(exp(y[1]+log(1/3))+exp(y[2]+log(1/3))+exp(y[3]+log(1/3)))<=0)
  JuMP.@NLconstraint(P1, log(exp(y[1]) + exp(y[2]) + exp(y[3])) <= log(3))
  #JuMP.@constraint(P1,y[1]>=log(0.1))
  #JuMP.@constraint(P1,y[2]>=log(0.1))
  #JuMP.@constraint(P1,y[3]>=log(0.1))
  JuMP.@NLobjective(P1, Min, 0.5 * (y[1])^2 + 0.5 * (y[2])^2)
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, y[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, 2 * y[1] + 2 * y[3] >= log(0.5))
  JuMP.@constraint(P2, 2 * y[2] + 2 * y[3] >= log(0.5))
  #JuMP.@constraint(P2,y[1]>=log(0.1))
  #JuMP.@constraint(P2,y[2]>=log(0.1))
  #JuMP.@constraint(P2,y[3]>=log(0.1))
  JuMP.@NLobjective(P2, Min, y[1] + y[2] + 2 * y[3])
  Player2 = MathProgNLPModel(P2)
  sol = x -> ((x[1] == x[2] && x[1] == (0.5 + 1 / sqrt(2))) || x == [1.0, 1.0, sqrt(2) / 2])
  #norm(xk-log.([1/2+1/sqrt(2),1/2+1/sqrt(2),1/(sqrt(2)*(1/2+1/sqrt(2)))]),Inf)
  #norm(xk-log.([1.,1.,sqrt(2)/2]),Inf)
  return GNEP(2, [2, 1], [Player1, Player2], sol = sol)
end

function Jason3()
  ul(i) = [0; 0][i]
  ux(i) = [Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1]^2 + x[2]^2 - 1 <= 0)
  #JuMP.@constraint(P1,x[1]>=0)
  #JuMP.@constraint(P1,x[2]>=0)
  #JuMP.@NLobjective(P1,Min,x[1]^2-x[1]*x[2]-x[2])
  JuMP.@NLobjective(P1, Min, x[1]^2 - x[1] * x[2] - x[1] + x[2]^2)
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1]^2 + x[2]^2 - 1 <= 0)
  #JuMP.@constraint(P2,x[1]>=0)
  #JuMP.@constraint(P2,x[2]>=0)
  #JuMP.@NLobjective(P2,Min,x[2]^2-x[1]*x[2]/2-2*x[2])
  JuMP.@NLobjective(P2, Min, x[2]^2 - x[1] * x[2] / 2 - 2 * x[2] + x[1]^2)
  Player2 = MathProgNLPModel(P2)
  sol = x -> (t = x[1]; x[2] == sqrt(1 - t^2) && t <= 4 / 5 && t >= 0)

  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function TangiCEx()
  P1 = JuMP.Model()
  JuMP.@variable(P1, x[1:4], start = 1.0)
  JuMP.@constraint(P1, x[1] + x[3] >= 0)
  JuMP.@constraint(P1, x[2] + x[4] >= 0)
  JuMP.@NLobjective(P1, Min, x[1] + 2 * x[2])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[1:4], start = 1.0)
  JuMP.@constraint(P2, x[1] + x[3] >= 0)
  JuMP.@constraint(P2, x[2] + x[4] >= 0)
  JuMP.@NLobjective(P2, Min, 2 * x[3] + 4 * x[4])
  Player2 = MathProgNLPModel(P2)
  sol = x -> true
  #0 is a solution / may be unbounded!!

  return GNEP(2, [2, 2], [Player1, Player2], sol = sol)
end

function GNEP1()
  ul(i) = [0; 0][i]
  ux(i) = [Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1]^2 + x[2]^2 - 1 <= 0)
  JuMP.@constraint(P1, x[1] + x[2] - 1 <= 0)
  JuMP.@NLobjective(P1, Min, 0.5 * (x[1] - 1)^2)
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + x[2] - 1 <= 0)
  JuMP.@NLobjective(P2, Min, 0.5 * (x[2] - 0.5)^2)
  Player2 = MathProgNLPModel(P2)
  sol = x -> (t = x[1]; x[2] == 1 - t && 0.5 <= t && t <= 1)
  #VE: (3/4,1/4)

  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function AbsVal()
  P1 = JuMP.Model()
  JuMP.@variable(P1, x[1:2], start = 1.0)
  JuMP.@constraint(P1, x[1] >= -1)
  JuMP.@constraint(P1, x[1] <= 1)
  JuMP.@NLobjective(P1, Min, 0.0)
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[1:2], start = 1.0)
  JuMP.@constraint(P2, x[1] >= -1)
  JuMP.@constraint(P2, x[1] <= 1)
  JuMP.@constraint(P2, x[2] >= -x[1])
  JuMP.@constraint(P2, x[2] <= x[1])
  JuMP.@NLobjective(P2, Min, -x[2])
  Player2 = MathProgNLPModel(P2)
  sol = x -> (t = x[1]; x[2] == abs(t) && -1 <= t && t <= 1)
  #VE: (3/4,1/4)

  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function A8()
  ul(i) = [0; 0; 0][i]
  ux(i) = [1; 1; 2][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + x[2] - 1 <= 0)
  JuMP.@constraint(P1, x[1] + x[2] - x[3] >= 0)
  #JuMP.@constraint(P1,x[1]>=0)
  #JuMP.@constraint(P1,x[2]>=0)
  #JuMP.@constraint(P1,x[3]>=0)
  #JuMP.@constraint(P1,x[3]<=2)
  JuMP.@NLobjective(P1, Min, -x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + x[2] - 1 <= 0)
  JuMP.@constraint(P2, x[1] + x[2] - x[3] >= 0)
  #JuMP.@constraint(P1,x[1]>=0)
  #JuMP.@constraint(P1,x[2]>=0)
  #JuMP.@constraint(P1,x[3]>=0)
  #JuMP.@constraint(P1,x[3]<=2)
  JuMP.@NLobjective(P2, Min, 0.5 * (x[2] - 0.5)^2)
  Player2 = MathProgNLPModel(P2)

  P3 = JuMP.Model()
  JuMP.@variable(P3, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  #JuMP.@constraint(P1,x[1]>=0)
  #JuMP.@constraint(P1,x[2]>=0)
  #JuMP.@constraint(P1,x[3]>=0)
  #JuMP.@constraint(P1,x[3]<=2)
  JuMP.@NLobjective(P3, Min, 0.5 * (x[3] - 1.5 * x[1])^2)
  Player3 = MathProgNLPModel(P3)
  sol = x -> (t = x[1]; x[2] == 1 - t && x[3] == 1.5 * t && 0.5 <= t && t <= 2 / 3)

  return GNEP(3, [1, 1, 1], [Player1, Player2, Player3], sol = sol)
end

function A8slack()
  ul(i) = [0; 0; 0; 0; 0; 0; 0][i]
  ux(i) = [1; 1; 2; Inf; Inf; Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:7], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + x[4] - 1 == -x[2])
  JuMP.@constraint(P1, x[1] + x[4] - x[7] == x[3])
  JuMP.@NLobjective(P1, Min, -x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:7], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + x[4] - 1 == -x[5])
  JuMP.@constraint(P2, x[1] + x[4] - x[7] == x[6])
  JuMP.@NLobjective(P2, Min, 0.5 * (x[2] - 0.5)^2)
  Player2 = MathProgNLPModel(P2)

  P3 = JuMP.Model()
  JuMP.@variable(P3, x[i = 1:7], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@NLobjective(P3, Min, 0.5 * (x[7] - 1.5 * x[1])^2)
  Player3 = MathProgNLPModel(P3)
  sol = x -> (t = x[1]; x[2] == 1 - t && x[3] == 1.5 * t && 0.5 <= t && t <= 2 / 3)

  return GNEP(3, [3, 3, 1], [Player1, Player2, Player3], sol = sol)
end

function Riverbasin()
  a1 = 0.01
  a2 = 0.05
  a3 = 0.01
  b = 0.01
  e1 = 2.9
  e2 = 2.88
  e3 = 2.85
  ul(i) = [0; 0; 0][i]
  ux(i) = [Inf; Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, 3.25 * x[1] + 1.25 * x[2] + 4.125 * x[3] <= 100)
  JuMP.@constraint(P1, 2.29115 * x[1] + 1.5625 * x[2] + 2.8125 * x[3] <= 100)
  JuMP.@NLobjective(P1, Min, (a1 * x[1] + b * (x[1] + x[2] + x[3]) - e1) * x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, 3.25 * x[1] + 1.25 * x[2] + 4.125 * x[3] <= 100)
  JuMP.@constraint(P2, 2.29115 * x[1] + 1.5625 * x[2] + 2.8125 * x[3] <= 100)
  JuMP.@NLobjective(P2, Min, (a2 * x[2] + b * (x[1] + x[2] + x[3]) - e2) * x[2])
  Player2 = MathProgNLPModel(P2)

  P3 = JuMP.Model()
  JuMP.@variable(P3, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P3, 3.25 * x[1] + 1.25 * x[2] + 4.125 * x[3] <= 100)
  JuMP.@constraint(P3, 2.29115 * x[1] + 1.5625 * x[2] + 2.8125 * x[3] <= 100)
  JuMP.@NLobjective(P3, Min, (a3 * x[3] + b * (x[1] + x[2] + x[3]) - e3) * x[3])
  Player3 = MathProgNLPModel(P3)
  sol = x -> (t = x[1]; x[2] == 1 - t && x[3] == 1.5 * t && 0.5 <= t && t <= 2 / 3)

  return GNEP(3, [1, 1, 1], [Player1, Player2, Player3], sol = sol)
end

function Harker()
  ul(i) = [0; 0][i]
  ux(i) = [10; 10][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + x[2] <= 15)
  JuMP.@NLobjective(P1, Min, x[1]^2 + 8 / 3 * x[1] * x[2] - 34 * x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + x[2] <= 15)
  JuMP.@NLobjective(P2, Min, x[2]^2 + 5 / 4 * x[1] * x[2] - 24.25 * x[2])
  Player2 = MathProgNLPModel(P2)

  sol = x -> (t = x[1]; (x[1] == 5 && x[2] == 9) || (x[2] == 15 - t && t <= 10 && t >= 9))
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function ContrHarker()
  ul(i) = [0; 0][i]
  ux(i) = [10; 10][i]

  theta = 0.9

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + 0.9 * x[2] <= 14.4)
  JuMP.@NLobjective(P1, Min, x[1]^2 + 8 / 3 * x[1] * x[2] - 34 * x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, 0.9 * x[1] + x[2] <= 14.1)
  JuMP.@NLobjective(P2, Min, x[2]^2 + 5 / 4 * x[1] * x[2] - 24.25 * x[2])
  Player2 = MathProgNLPModel(P2)

  sol = x -> (t = x[1]; (x[1] == 5 && x[2] == 9) || (x[2] == 15 - t && t <= 10 && t >= 9))
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function Nabutani1()
  ul(i) = [0; 0; 0][i]
  ux(i) = [Inf; Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + 2 * x[2] - x[3] <= 14)
  JuMP.@constraint(P1, 3 * x[1] + 2 * x[2] + x[3] <= 30)
  JuMP.@NLobjective(
    P1,
    Min,
    x[1]^2 + x[1] * x[2] + x[2]^2 + (x[1] + x[2]) * x[3] - 25 * x[1] - 38 * x[2]
  )
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + 2 * x[2] - x[3] <= 14)
  JuMP.@constraint(P2, 3 * x[1] + 2 * x[2] + x[3] <= 30)
  JuMP.@NLobjective(P2, Min, x[3]^2 + (x[1] + x[2]) * x[3] - 25 * x[3])
  Player2 = MathProgNLPModel(P2)

  sol = x -> (t = x[1]; x[2] == 11 - t && x[3] == 8 - t && t >= 0 && t <= 2)
  return GNEP(2, [2, 1], [Player1, Player2], sol = sol)
end

function Ex1()
  ul(i) = [0; 0][i]
  ux(i) = [Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + x[2] <= 1)
  JuMP.@NLobjective(P1, Min, x[1]^2 - x[1] * x[2] - x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + x[2] <= 1)
  JuMP.@NLobjective(P2, Min, x[2]^2 - 1 / 2 * x[1] * x[2] - 2 * x[2])
  Player2 = MathProgNLPModel(P2)

  sol = x -> (t = x[1]; x[2] == 1 - t && t <= 2 / 3 && t >= 0)
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function Ex1bis()
  ul(i) = [0; 0][i]
  ux(i) = [Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1]^2 + x[2]^2 <= 1)
  JuMP.@NLobjective(P1, Min, x[1]^2 - x[1] * x[2] - x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1]^2 + x[2]^2 <= 1)
  JuMP.@NLobjective(P2, Min, x[2]^2 - 1 / 2 * x[1] * x[2] - 2 * x[2])
  Player2 = MathProgNLPModel(P2)

  sol = x -> (t = x[1]; x[2] == sqrt(1 - t^2) && t <= 4 / 5 && t >= 0)
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function ContreExDS()
  ul(i) = [-Inf; -Inf][i]
  ux(i) = [Inf; Inf][i]
  L = 2

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] <= 1)
  JuMP.@constraint(P1, -x[1] + L * x[2] <= 0)
  JuMP.@NLobjective(P1, Min, x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] - x[2] == 0)
  JuMP.@NLobjective(P2, Min, (x[1] - x[2])^2)
  Player2 = MathProgNLPModel(P2)

  sol = x -> (t = x[1]; x[2] == t && t == L * x[2]) #no solutions if L>1, (1,1) if L=1, any solution if L<1
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function Ex4()
  ul(i) = [0.0; -Inf; 0.0; -Inf][i]
  ux(i) = [Inf; Inf; Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:4], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + x[3] <= 1)
  JuMP.@constraint(P1, x[1] - 2 * x[3] - x[2] <= 0)
  JuMP.@constraint(P1, -x[1] - x[3] - x[2] <= 0)
  JuMP.@NLobjective(P1, Min, x[2])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:4], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + x[3] <= 1)
  JuMP.@constraint(P2, -x[1] + x[3] - x[4] <= 0)
  JuMP.@constraint(P2, x[1] - x[3] - x[4] <= 0)
  JuMP.@NLobjective(P2, Min, x[4])
  Player2 = MathProgNLPModel(P2)

  sol = x -> (true)#((x[3]=1-x[1]) || ((x[3]==x[1]) && (x[1]<=0.5)))
  return GNEP(2, [2, 2], [Player1, Player2], sol = sol)
end

function Ex5()
  ul(i) = [-Inf; -Inf; 0][i]
  ux(i) = [Inf; Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1] + x[2] + x[3] <= 1)
  JuMP.@NLobjective(P1, Min, 0.5 * (x[1] - 1)^2 - x[1] * x[2])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P2, x[1] + x[2] + x[3] <= 1)
  JuMP.@NLobjective(P2, Min, 0.5 * (x[2] - 1)^2 + x[1] * x[2])
  Player2 = MathProgNLPModel(P2)

  P3 = JuMP.Model()
  JuMP.@variable(P3, x[i = 1:3], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P3, -x[1] - x[2] + x[3] <= 0)
  JuMP.@NLobjective(P3, Min, 0.5 * (x[3] - 1)^2)
  Player3 = MathProgNLPModel(P3)

  sol = x -> (true)

  return GNEP(3, [1, 1, 1], [Player1, Player2, Player3], sol = sol)
end

function Ex6()
  ul(i) = [-Inf; -Inf][i]
  ux(i) = [Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, 1 / 6 * x[1]^2 + x[2] - 5 / 2 <= 0)
  JuMP.@NLobjective(P1, Min, 0.5 * x[1]^2 + 32 / 5 * x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@NLobjective(P2, Min, 0.5 * x[2]^2 + x[1] * x[2] - 4 / 5 * x[2])
  Player2 = MathProgNLPModel(P2)

  sol = x -> (true)
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function Ex7()
  ul(i) = [-Inf; -Inf][i]
  ux(i) = [Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1]^2 + x[2] <= 1)
  JuMP.@NLobjective(P1, Min, x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@NLobjective(P2, Min, 0.5 * x[2]^2)
  Player2 = MathProgNLPModel(P2)

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

function Ex8()
  ul(i) = [-Inf; -Inf][i]
  ux(i) = [Inf; Inf][i]

  P1 = JuMP.Model()
  JuMP.@variable(P1, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@constraint(P1, x[1]^2 + x[2] <= 0)
  JuMP.@NLobjective(P1, Min, 0.5 * x[1]^2 - 2 * x[1])
  Player1 = MathProgNLPModel(P1)

  P2 = JuMP.Model()
  JuMP.@variable(P2, x[i = 1:2], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
  JuMP.@NLobjective(P2, Min, 0.5 * x[2]^2 + (2 - x[1]^2) * x[2])
  Player2 = MathProgNLPModel(P2)

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(2, [1, 1], [Player1, Player2], sol = sol)
end

"""
Arrow and Debreu model of a competitive economy proposed in:
Kenneth J. Arrow and Gerard Debreu. Existence of an equilibrium for a competitive
economy. Econometrica: Journal of the Econometric Society, pages 265–290, 1954.

Also known as problem A10 in:
Francisco Facchinei and Christian Kanzow. Penalty methods for the solution of general-
ized Nash equilibrium problems (with complete test problems). Sapienza University of
Rome, 2009.

This problem is specialized in several variants: a), b), c), d) and e)
"""
function A10a()
  include("A10a.jl")
  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, F + C + 1)
  ul(i) = zeros((F + C + 1) * P)[i]
  ux(i) = vcat(sqrt(21) * ones(P * F), 6 * ones(P * C), 1.1 * ones(P))[i]

  #Firms problems
  for j = 1:F
    firm_j = JuMP.Model()
    JuMP.@variable(
      firm_j,
      x[i = 1:((F + C + 1) * P)],
      upperbound = ux(i),
      lowerbound = ul(i),
      start = 1.0
    )
    #JuMP.@NLobjective(firm_j,Min,x[(F+C)*P+1:n]'*x[P*(j-1)+1:P*(j-1)+P])
    JuMP.@NLobjective(
      firm_j,
      Min,
      -x[22] * x[P * (j - 1) + 1] - x[23] * x[P * (j - 1) + 2] - x[24] * x[P * (j - 1) + 3]
    )
    JuMP.@constraint(firm_j, sum(x[(P * (j - 1) + 1):(P * (j - 1) + P)] .^ 2) <= 10 * j)
    arr_player[j] = MathProgNLPModel(firm_j)
  end
  #Consumers problems
  for i = 1:C
    cons_i = JuMP.Model()
    JuMP.@variable(
      cons_i,
      x[i = 1:((F + C + 1) * P)],
      upperbound = ux(i),
      lowerbound = ul(i),
      start = 1.0
    )
    #error max not min
    #JuMP.@NLobjective(cons_i,Min, -0.5*(x[1]*Q[i,1,1]*x[1]+x[1]*Q[i,1,2]*x[2]+x[1]*Q[i,1,1]*x[3] + x[2]*Q[i,2,1]*x[1]+x[2]*Q[i,2,2]*x[2]+x[2]*Q[i,2,1]*x[3] + x[3]*Q[i,3,1]*x[1]+x[3]*Q[i,3,2]*x[2]+x[3]*Q[i,3,1]*x[3]) + x[P*(F+i-1)+1]*b[i,1,1] + x[P*(F+i-1)+2]*b[i,2,1] +x[P*(F+i-1)+3]*b[i,3,1])
    JuMP.@NLobjective(
      cons_i,
      Min,
      0.5 * (
        x[1] * Q[i, 1, 1] * x[1] +
        x[1] * Q[i, 1, 2] * x[2] +
        x[1] * Q[i, 1, 1] * x[3] +
        x[2] * Q[i, 2, 1] * x[1] +
        x[2] * Q[i, 2, 2] * x[2] +
        x[2] * Q[i, 2, 1] * x[3] +
        x[3] * Q[i, 3, 1] * x[1] +
        x[3] * Q[i, 3, 2] * x[2] +
        x[3] * Q[i, 3, 1] * x[3]
      ) - x[P * (F + i - 1) + 1] * b[i, 1, 1] - x[P * (F + i - 1) + 2] * b[i, 2, 1] -
      x[P * (F + i - 1) + 3] * b[i, 3, 1]
    )
    JuMP.@constraint(
      cons_i,
      x[((F + C) * P + 1):((F + C + 1) * P)]' * x[(P * (F + i - 1) + 1):(P * (F + i - 1) + P)] -
      x[((F + C) * P + 1):((F + C + 1) * P)]' * xi[i, :, 1] <= 0
    )
    arr_player[F + i] = MathProgNLPModel(cons_i)
  end
  #Market player problem
  mark = JuMP.Model()
  JuMP.@variable(
    mark,
    x[i = 1:((F + C + 1) * P)],
    upperbound = ux(i),
    lowerbound = ul(i),
    start = 1.0
  )
  JuMP.@NLobjective(
    mark,
    Min,
    -x[22] * (
      x[7] + x[10] + x[13] + x[16] + x[19] - x[1] - x[4] - xi[1, 1, 1] - xi[2, 1, 1] - xi[3, 1, 1] -
      xi[4, 1, 1] - xi[5, 1, 1]
    ) -
    x[23] * (
      x[8] + x[11] + x[14] + x[17] + x[20] - x[2] - x[5] - xi[1, 2, 1] - xi[2, 2, 1] - xi[3, 2, 1] -
      xi[4, 2, 1] - xi[5, 2, 1]
    ) -
    x[24] * (
      x[9] + x[12] + x[15] + x[18] + x[21] - x[3] - x[5] - xi[1, 3, 1] - xi[2, 3, 1] - xi[3, 3, 1] -
      xi[4, 3, 1] - xi[5, 3, 1]
    )
  )
  JuMP.@constraint(mark, sum(x[((F + C) * P + 1):((F + C + 1) * P)]) == 1)
  arr_player[F + C + 1] = MathProgNLPModel(mark)

  sol = x -> (true)

  #initial point in the paper:
  x0 = zeros(24)
  x0[22:24] .= 1 / 3
  return GNEPmod.GNEP(F + C + 1, P * ones(Int64, F + C + 1), arr_player, sol = sol, x0 = x0)
end

"""
Inspired by an internet switching model introduced by

Alex Kesselman, Stefano Leonardi, and Vincenzo Bonifaci. Game-theoretic analysis of
internet switching with selfish users. In Xiaotie Deng and Yinyu Ye, editors, Internet and
Network Economics, pages 236–245, Berlin, Heidelberg, 2005. Springer Berlin Heidelberg.

and also analyzed as problem A14 in
Francisco Facchinei and Christian Kanzow. Penalty methods for the solution of general-
ized Nash equilibrium problems (with complete test problems). Sapienza University of
Rome, 2009.

This is my* specialization.
"""
function TwoGroupsSwitch(B1, B2)
  N1 = 2
  N2 = 3
  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N1 + N2)
  ul(i) = 0.01 * ones(Int64, N1 + N2)[i]
  ux(i) = Inf * ones(Int64, N1 + N2)[i]

  for k = 1:N1
    P1 = JuMP.Model()
    JuMP.@variable(P1, x[i = 1:(N1 + N2)], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, sum(x[1:(N1 + N2)]) <= B1)
    JuMP.@NLobjective(P1, Min, -x[k] / (x[1] + x[2] + x[3] + x[4] + x[5]))
    arr_player[k] = MathProgNLPModel(P1)
  end

  for j = (N1 + 1):(N1 + N2)
    P2 = JuMP.Model()
    JuMP.@variable(P2, x[i = 1:(N1 + N2)], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P2, sum(x[1:(N1 + N2)] .^ 2) <= B2)
    JuMP.@NLobjective(P2, Min, -x[j] / (x[1] + x[2] + x[3] + x[4] + x[5]))
    arr_player[j] = MathProgNLPModel(P2)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N1 + N2, ones(Int64, N1 + N2), arr_player, sol = sol)
end

using LinearAlgebra
"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.
"""
function EPC()
  include("EPC.jl")

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ub = 50 * ones(Float64, N * (N + 1))
  ub[union(1, 5, 9)] = p

  ul(i) = zeros(Float64, N * (N + 1))[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * (N + 1) #index of e_k
    j = mod(k, N) + 1
    l = mod(k + 1, N) + 1
    JuMP.@variable(P1, x[i = 1:(N * (N + 1))], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, dot(x[nu:(N + nu)], vcat(1.0, -gam)) <= E[k])
    JuMP.@constraint(P1, -dot(x[union(1, 2, 6, 10)], vcat(1.0, -gam[1] * ones(N))) <= 0.0)
    JuMP.@constraint(P1, -dot(x[union(5, 3, 7, 11)], vcat(1.0, -gam[2] * ones(N))) <= 0.0)
    JuMP.@constraint(P1, -dot(x[union(9, 4, 8, 12)], vcat(1.0, -gam[3] * ones(N))) <= 0.0)
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (
        x[nu + k]^2 +
        x[nu + j]^2 +
        x[nu + l]^2 +
        2 * x[nu + j] * x[1 + (j - 1) * (N + 1) + j] +
        2 * x[nu + l] * x[1 + (l - 1) * (N + 1) + l]
      ) +
      lambda[1] * (x[1] - gam[1] * (x[2] + x[6] + x[10])) +
      lambda[2] * (x[5] - gam[2] * (x[3] + x[7] + x[11])) +
      lambda[3] * (x[9] - gam[3] * (x[4] + x[8] + x[12]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, (N + 1) * ones(Int64, N), arr_player, sol = sol)
end

"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.

Variant of EPC() studied in
Cojocaru, M. G., Wild, E., & Small, A. (2018).
On describing the solution sets of generalized Nash games with shared constraints.
Optimization and Engineering, 19(4), 845-870.
"""
function EPCshared()
  include("EPC.jl")

  ub = 50 * ones(Float64, N * (N + 1))
  ub[union(1, 5, 9)] = p

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = zeros(Float64, N * (N + 1))[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * (N + 1) #index of e_k
    j = mod(k, N) + 1
    l = mod(k + 1, N) + 1
    JuMP.@variable(P1, x[i = 1:(N * (N + 1))], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, dot(x[nu:(N + nu)], vcat(1.0, -gam)) <= E[k])
    JuMP.@constraint(
      P1,
      -dot(x[union(1, 2, 6, 10)], vcat(1.0, -gam[1] * ones(N))) -
      dot(x[union(5, 3, 7, 11)], vcat(1.0, -gam[2] * ones(N))) -
      dot(x[union(9, 4, 8, 12)], vcat(1.0, -gam[3] * ones(N))) <= 0.0
    )
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (
        x[nu + k]^2 +
        x[nu + j]^2 +
        x[nu + l]^2 +
        2 * x[nu + j] * x[1 + (j - 1) * (N + 1) + j] +
        2 * x[nu + l] * x[1 + (l - 1) * (N + 1) + l]
      ) +
      lambda[1] * (x[1] - gam[1] * (x[2] + x[6] + x[10])) +
      lambda[2] * (x[5] - gam[2] * (x[3] + x[7] + x[11])) +
      lambda[3] * (x[9] - gam[3] * (x[4] + x[8] + x[12]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, (N + 1) * ones(Int64, N), arr_player, sol = sol)
end

"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.

Variant of EPC() to find the emission of each country
"""
function regularEPC()
  include("EPC.jl")
  ub = 100 * ones(Float64, 6)

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = zeros(Int64, 2 * N)[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * 2 #index of e_k
    JuMP.@variable(P1, x[i = 1:(2 * N)], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, x[nu] - gam[k] * x[nu + 1] <= E[k])
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (x[nu + 1]^2) +
      lambda[1] * (x[1] - gam[1] * (x[2])) +
      lambda[2] * (x[3] - gam[2] * (x[4])) +
      lambda[3] * (x[5] - gam[3] * (x[6]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, 2 * ones(Int64, N), arr_player, sol = sol)
end

"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.

Variant of EPC() studied in *** new ***.
"""
function EPCnonshared()
  include("EPC.jl")

  ub = 50 * ones(Float64, N * (N + 1))
  ub[union(1, 5, 9)] = p

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = zeros(Float64, N * (N + 1))[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * (N + 1) #index of e_k
    j = mod(k, N) + 1
    l = mod(k + 1, N) + 1
    JuMP.@variable(P1, x[i = 1:(N * (N + 1))], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, dot(x[nu:(N + nu)], vcat(1.0, -gam)) <= E[k])
    JuMP.@constraint(
      P1,
      -dot(x[union(1, 2, 6, 10)], vcat(1.0, -gam[1] * ones(N))) -
      dot(x[union(5, 3, 7, 11)], vcat(1.0, -gam[2] * ones(N))) -
      dot(x[union(9, 4, 8, 12)], vcat(1.0, -gam[3] * ones(N))) <= 0.0
    )
    if (k == 1) || (k == 2)
      JuMP.@NLconstraint(
        P1,
        -(x[2] + x[6]) / (x[2] + x[3] + x[4] + x[6]) -
        (x[3] + x[7]) / (x[6] + x[7] + x[8] + x[3]) <= -gap
      )
    end
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (
        x[nu + k]^2 +
        x[nu + j]^2 +
        x[nu + l]^2 +
        2 * x[nu + j] * x[1 + (j - 1) * (N + 1) + j] +
        2 * x[nu + l] * x[1 + (l - 1) * (N + 1) + l]
      ) +
      lambda[1] * (x[1] - gam[1] * (x[2] + x[6] + x[10])) +
      lambda[2] * (x[5] - gam[2] * (x[3] + x[7] + x[11])) +
      lambda[3] * (x[9] - gam[3] * (x[4] + x[8] + x[12]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, (N + 1) * ones(Int64, N), arr_player, sol = sol)
end

"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.

Variant of EPC() studied in
Cojocaru, M. G., Wild, E., & Small, A. (2018).
On describing the solution sets of generalized Nash games with shared constraints.
Optimization and Engineering, 19(4), 845-870.

N=5
"""
function EPC5shared()
  include("EPC5.jl")

  ub = 50 * ones(Float64, N * (N + 1))
  ub[union(1, 7, 13, 19, 25)] = p

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = zeros(Float64, N * (N + 1))[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * (N + 1) #index of e_k
    j = mod(k, N) + 1
    l = mod(k + 1, N) + 1
    o = mod(k + 2, N) + 1
    m = mod(k + 3, N) + 1
    JuMP.@variable(P1, x[i = 1:(N * (N + 1))], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, dot(x[nu:(N + nu)], vcat(1.0, -gam)) <= E[k])
    JuMP.@constraint(
      P1,
      -dot(x[union(1, 2, 8, 14, 20, 26)], vcat(1.0, -gam[1] * ones(N))) -
      dot(x[union(7, 3, 9, 15, 21, 27)], vcat(1.0, -gam[2] * ones(N))) -
      dot(x[union(13, 4, 10, 16, 22, 28)], vcat(1.0, -gam[3] * ones(N))) -
      dot(x[union(19, 5, 11, 17, 23, 29)], vcat(1.0, -gam[4] * ones(N))) -
      dot(x[union(25, 6, 12, 18, 24, 30)], vcat(1.0, -gam[5] * ones(N))) <= 0.0
    )
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (
        x[nu + k]^2 +
        x[nu + j]^2 +
        x[nu + l]^2 +
        x[nu + o]^2 +
        x[nu + m]^2 +
        2 * x[nu + j] * x[1 + (j - 1) * (N + 1) + j] +
        2 * x[nu + l] * x[1 + (l - 1) * (N + 1) + l] +
        2 * x[nu + o] * x[1 + (o - 1) * (N + 1) + o] +
        2 * x[nu + m] * x[1 + (m - 1) * (N + 1) + m]
      ) +
      lambda[1] * (x[1] - gam[1] * (x[2] + x[8] + x[14] + x[20] + x[26])) +
      lambda[2] * (x[7] - gam[2] * (x[3] + x[9] + x[15] + x[21] + x[27])) +
      lambda[3] * (x[13] - gam[3] * (x[4] + x[10] + x[16] + x[22] + x[28])) +
      lambda[4] * (x[19] - gam[4] * (x[5] + x[11] + x[17] + x[23] + x[29])) +
      lambda[5] * (x[25] - gam[5] * (x[6] + x[12] + x[18] + x[24] + x[30]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, (N + 1) * ones(Int64, N), arr_player, sol = sol)
end

"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.

Variant of EPC() to find the emission of each country
"""
function regularEPC5()
  include("EPC5.jl")
  ub = 100 * ones(Float64, 2 * N)

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = zeros(Int64, 2 * N)[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * 2 #index of e_k
    JuMP.@variable(P1, x[i = 1:(2 * N)], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, x[nu] - gam[k] * x[nu + 1] <= E[k])
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (x[nu + 1]^2) +
      lambda[1] * (x[1] - gam[1] * (x[2])) +
      lambda[2] * (x[3] - gam[2] * (x[4])) +
      lambda[3] * (x[5] - gam[3] * (x[6])) +
      lambda[4] * (x[7] - gam[4] * (x[8])) +
      lambda[5] * (x[9] - gam[5] * (x[10]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, 2 * ones(Int64, N), arr_player, sol = sol)
end

"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.

Variant of EPC() studied in *** new ***.
"""
function EPC5nonshared()
  include("EPC5.jl")

  ub = 50 * ones(Float64, N * (N + 1))
  ub[union(1, 7, 13, 19, 25)] = p

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = zeros(Float64, N * (N + 1))[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * (N + 1) #index of e_k
    j = mod(k, N) + 1
    l = mod(k + 1, N) + 1
    o = mod(k + 2, N) + 1
    m = mod(k + 3, N) + 1
    JuMP.@variable(P1, x[i = 1:(N * (N + 1))], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, dot(x[nu:(N + nu)], vcat(1.0, -gam)) <= E[k])
    JuMP.@constraint(
      P1,
      -dot(x[union(1, 2, 8, 14, 20, 26)], vcat(1.0, -gam[1] * ones(N))) -
      dot(x[union(7, 3, 9, 15, 21, 27)], vcat(1.0, -gam[2] * ones(N))) -
      dot(x[union(13, 4, 10, 16, 22, 28)], vcat(1.0, -gam[3] * ones(N))) -
      dot(x[union(19, 5, 11, 17, 23, 29)], vcat(1.0, -gam[4] * ones(N))) -
      dot(x[union(25, 6, 12, 18, 24, 30)], vcat(1.0, -gam[5] * ones(N))) <= 0.0
    )
    if (k == 3) || (k == 4)
      JuMP.@NLconstraint(
        P1,
        -(x[16] * gam[4]) / (x[14] + x[15] + x[16] + x[17] + x[18]) -
        (x[23] * gam[3]) / (x[20] + x[21] + x[22] + x[23] + x[24]) <= -gap
      )
    end
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (
        x[nu + k]^2 +
        x[nu + j]^2 +
        x[nu + l]^2 +
        x[nu + o]^2 +
        x[nu + m]^2 +
        2 * x[nu + j] * x[1 + (j - 1) * (N + 1) + j] +
        2 * x[nu + l] * x[1 + (l - 1) * (N + 1) + l] +
        2 * x[nu + o] * x[1 + (o - 1) * (N + 1) + o] +
        2 * x[nu + m] * x[1 + (m - 1) * (N + 1) + m]
      ) +
      lambda[1] * (x[1] - gam[1] * (x[2] + x[8] + x[14] + x[20] + x[26])) +
      lambda[2] * (x[7] - gam[2] * (x[3] + x[9] + x[15] + x[21] + x[27])) +
      lambda[3] * (x[13] - gam[3] * (x[4] + x[10] + x[16] + x[22] + x[28])) +
      lambda[4] * (x[19] - gam[4] * (x[5] + x[11] + x[17] + x[23] + x[29])) +
      lambda[5] * (x[25] - gam[5] * (x[6] + x[12] + x[18] + x[24] + x[30]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, (N + 1) * ones(Int64, N), arr_player, sol = sol)
end

"""
Environmental protocol control proposed in
Breton, M., Zaccour, G., & Zahaf, M. (2006).
A game-theoretic formulation of joint implementation of environmental projects.
European Journal of Operational Research, 168(1), 221-239.

Variant of EPC() studied in *** new ***.
"""
function EPC5nonsharedgreen()
  include("EPC5.jl")

  ub = 50 * ones(Float64, N * (N + 1))
  ub[union(1, 7, 13, 19, 25)] = p

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = zeros(Float64, N * (N + 1))[i]
  ux(i) = ub[i]

  for k = 1:N
    P1 = JuMP.Model()
    nu = 1 + (k - 1) * (N + 1) #index of e_k
    j = mod(k, N) + 1
    l = mod(k + 1, N) + 1
    o = mod(k + 2, N) + 1
    m = mod(k + 3, N) + 1
    JuMP.@variable(P1, x[i = 1:(N * (N + 1))], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(P1, dot(x[nu:(N + nu)], vcat(1.0, -gam)) <= E[k])
    JuMP.@constraint(
      P1,
      -dot(x[union(1, 2, 8, 14, 20, 26)], vcat(1.0, -gam[1] * ones(N))) -
      dot(x[union(7, 3, 9, 15, 21, 27)], vcat(1.0, -gam[2] * ones(N))) -
      dot(x[union(13, 4, 10, 16, 22, 28)], vcat(1.0, -gam[3] * ones(N))) -
      dot(x[union(19, 5, 11, 17, 23, 29)], vcat(1.0, -gam[4] * ones(N))) -
      dot(x[union(25, 6, 12, 18, 24, 30)], vcat(1.0, -gam[5] * ones(N))) <= 0.0
    )
    if (k == 1)
      #JuMP.@NLconstraint(P1,x[3]-1/(1e-3+x[7]-97.5)^2<=0.0)
      #JuMP.@NLconstraint(P1,x[4]-1/(1e-3+x[13]-97.5)^2<=0.0)
      #JuMP.@NLconstraint(P1,x[5]-1/(1e-3+x[19]-97.5)^2<=0.0)
      #JuMP.@NLconstraint(P1,x[6]-1/(1e-3+x[25]-97.5)^2<=0.0)
      JuMP.@NLconstraint(P1, x[3] * (x[7] - 97.5)^2 - 1 <= 0.0)
      JuMP.@NLconstraint(P1, x[4] * (x[13] - 97.5)^2 - 1 <= 0.0)
      JuMP.@NLconstraint(P1, x[5] * (x[19] - 97.5)^2 - 1 <= 0.0)
      JuMP.@NLconstraint(P1, x[6] * (x[25] - 97.5)^2 - 1 <= 0.0)
    end
    JuMP.@NLobjective(
      P1,
      Min,
      -x[nu] * (p[k] - 0.5 * x[nu]) +
      0.5 * (
        x[nu + k]^2 +
        x[nu + j]^2 +
        x[nu + l]^2 +
        x[nu + o]^2 +
        x[nu + m]^2 +
        2 * x[nu + j] * x[1 + (j - 1) * (N + 1) + j] +
        2 * x[nu + l] * x[1 + (l - 1) * (N + 1) + l] +
        2 * x[nu + o] * x[1 + (o - 1) * (N + 1) + o] +
        2 * x[nu + m] * x[1 + (m - 1) * (N + 1) + m]
      ) +
      lambda[1] * (x[1] - gam[1] * (x[2] + x[8] + x[14] + x[20] + x[26])) +
      lambda[2] * (x[7] - gam[2] * (x[3] + x[9] + x[15] + x[21] + x[27])) +
      lambda[3] * (x[13] - gam[3] * (x[4] + x[10] + x[16] + x[22] + x[28])) +
      lambda[4] * (x[19] - gam[4] * (x[5] + x[11] + x[17] + x[23] + x[29])) +
      lambda[5] * (x[25] - gam[5] * (x[6] + x[12] + x[18] + x[24] + x[30]))
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, (N + 1) * ones(Int64, N), arr_player, sol = sol)
end

"""
A scalable game
"""
function NeighborsGame(N::Int64)
  if N <= 0
    throw("Error: N must be positive.")
  end

  lb = zeros(Float64, N)
  ub = 10 * ones(Float64, N)

  arr_player = Array{NLPModelsJuMP.MathProgNLPModel, 1}(undef, N)

  ul(i) = lb[i]
  ux(i) = ub[i]

  C = 0.5 * rand(N, 4)

  for k = 1:N
    km2 = mod(k - 3, N) + 1
    km1 = mod(k - 2, N) + 1
    kp1 = mod(k, N) + 1
    kp2 = mod(k + 1, N) + 1

    P1 = JuMP.Model()
    JuMP.@variable(P1, x[i = 1:N], upperbound = ux(i), lowerbound = ul(i), start = 1.0)
    JuMP.@constraint(
      P1,
      x[k]^2 + x[km2]^2 + x[km1]^2 + x[kp1]^2 + x[kp2]^2 - C[k, 1] * x[k] * x[km2] -
      C[k, 2] * x[k] * x[km1] - C[k, 3] * x[k] * x[kp1] - C[k, 4] * x[k] * x[kp2] <= 0
    )
    JuMP.@NLobjective(
      P1,
      Min,
      x[k]^2 + x[km2]^2 + x[km1]^2 + x[kp1]^2 + x[kp2]^2 - x[k] * x[kp1] * x[kp2] -
      x[k] * x[km1] * x[km2]
    )
    arr_player[k] = MathProgNLPModel(P1)
  end

  sol = x -> ((x[1] == -1.0) && (x[2] == 0.0))
  return GNEP(N, ones(Int64, N), arr_player, sol = sol)
end


text:
# module Bukdu.Plug

struct Auth <: AbstractPlug
end

function plug(::Type{Auth}, conn::Conn)
    # TODO
    @error Auth @__FILE__
    conn.request.response.status = 401 # 401 Unauthorized
    conn.halted = true
end

# module Bukdu.Plug

ext: jl

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
using Turing.VarReplay

mutable struct Trace
  task  ::  Task
  vi    ::  VarInfo
  spl   ::  Union{Nothing, Sampler}
  Trace() = (res = new(); res.vi = VarInfo(); res.spl = nothing; res)
end

# NOTE: this function is called by `forkr`
function Trace(f)
  res = Trace();
  # Task(()->f());
  res.task = Task( () -> begin res=f(); produce(Val{:done}); res; end )
  if isa(res.task.storage, Nothing)
    res.task.storage = IdDict()
  end
  res.task.storage[:turing_trace] = res # create a backward reference in task_local_storage
  res
end

function Trace(f, spl::Sampler, vi :: VarInfo)
  res = Trace();
  res.spl = spl
  # Task(()->f());
  res.vi = deepcopy(vi)
  res.vi.num_produce = 0
  res.task = Task( () -> begin vi_new=f(vi, spl); produce(Val{:done}); vi_new; end )
  if isa(res.task.storage, Nothing)
    res.task.storage = IdDict()
  end
  res.task.storage[:turing_trace] = res # create a backward reference in task_local_storage
  res
end

# step to the next observe statement, return log likelihood
Turing.consume(t::Trace) = (t.vi.num_produce += 1; consume(t.task))

# Task copying version of fork for Trace.
function fork(trace :: Trace, is_ref :: Bool = false)
  newtrace = typeof(trace)()
  newtrace.task = Base.copy(trace.task)
  newtrace.spl = trace.spl

  newtrace.vi = deepcopy(trace.vi)
  if is_ref
    set_retained_vns_del_by_spl!(newtrace.vi, newtrace.spl)
  end

  newtrace.task.storage[:turing_trace] = newtrace
  newtrace
end

# PG requires keeping all randomness for the reference particle
# Create new task and copy randomness
function forkr(trace :: Trace)
  newtrace = Trace(trace.task.code)
  newtrace.spl = trace.spl

  newtrace.vi = deepcopy(trace.vi)
  newtrace.vi.num_produce = 0

  newtrace
end

current_trace() = current_task().storage[:turing_trace]


text:
# tests for the synchronous task runner
import Arbiter.Sync: run_tasks
import Arbiter.ArbiterTasks: ArbiterTask
import Arbiter.ArbiterGraphs: NodeSet, ImmutableNodeSet

facts("empty") do
    # solve no tasks

    results = run_tasks(())

    @fact results.completed --> ImmutableNodeSet()
    @fact results.failed --> ImmutableNodeSet()
end

facts("no dependencies") do
    # run dependency-less tasks

    executed_tasks = NodeSet()

    """
    Make a task
    """
    function make_task(name, dependencies=(); succeed=true)
        ArbiterTask(name, () -> (push!(executed_tasks, name); succeed), dependencies)
    end

    results = run_tasks(
        (
            make_task(:foo),
            make_task(:bar),
            make_task(:baz),
            make_task(:fail; succeed=false)
        )
    )

    @fact executed_tasks --> NodeSet((:foo, :bar, :baz, :fail))
    @fact results.completed --> ImmutableNodeSet((:foo, :bar, :baz))
    @fact results.failed --> ImmutableNodeSet((:fail,))
end


text:
struct SigMFDataTypeError <: Exception
    s::String
end

Base.showerror(io::IO, e::SigMFDataTypeError) = print(io, "Could not parse datatype format: ", e.s)


text:
module ExtractPSF

using OffsetArrays, Optim, Images, Statistics

import Base: size, getindex, setindex!, axes, parentindices

include("find.jl")
include("optics.jl")
include("fit_gaussian.jl")
include("psf.jl")
include("synthetic_psfs.jl")

export find_beads,
        Psf, gaussian_lightsheet_psf,
        PsfFit, psf, shift, quality, scale

end # module

max_forks_repo_name: HolyLab/ExtractPSF


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function unit_conversions_test(fname)
    pagewidth, pageheight = 2000, 2000
    Drawing(pagewidth, pageheight, fname)
    origin()
    background("ivory")
    setline(0.5)
    tiles = Tiler(pagewidth, pageheight, 10, 4, margin=10)
    p = rand()
    for (pos, n) in tiles
        gsave()
        sethue("blue")
        translate(pos)
        if n % 2 == 0
            rect(O, 2inch, 1inch, :stroke)
            circle(O, 3inch/10, :stroke)

            rect(O + Point(-5mm, 3mm), 2cm, 1cm, :stroke)
            circle(O, 1cm, :stroke)

            rect(O, p*cm, p*cm, :stroke)
            circle(O, p*cm, :stroke)

        else
            sethue("darkorange")
            rect(O, 2*72, 72, :stroke)
            circle(O, 216/10, :stroke)

            rect(O + Point(-5 * 2.83464566929, 3 * 2.83464566929), 2 * 28.3464566929, 28.3464566929, :stroke)
            circle(O, 28.3464566929, :stroke)

            rect(O, p * 28.3464566929, p * 28.3464566929, :stroke)
            circle(O, p * 28.3464566929, :stroke)
        end
        grestore()
    end
    @test finish() == true
end

@test isapprox(1inch, 72)
@test isapprox(72, 1inch)
@test isapprox(2.54cm, 1inch)
@test isapprox(72, 2.54cm)

fname = "unit-conversions-test.pdf"
unit_conversions_test(fname)
println("...finished test: output in $(fname)")


text:
export iteratedTikhonov

"""
    mc,DC = iteratedTikhonov(mc,pInv::InverseParam,pMis,nAlpha,alphaFac,
                            targetMisfit;indCredit=[],dumpResults::Function=dummy)

    Perform (Projected) Gauss-NewtonCG using iterated Tikhonov procedure to decrease
    regularization parameter and update reference model after fixed number of GN iterations
    set in pInv.

    Input:
           mc::Vector         - Initial guess for model
           pInv::InverseParam - parameter for inversion
           pMis               - misfit terms
           nAlpha             - maximum number of allowed regularization parameter (alpha) values
           alphaFac           - alpha decrease factor. alpha_(i+1) = alpha_i/alphaFac
           targetMisfit       - Termination criterion. iteratedTikhonov will exit
                                 -when data misfit < targetMisfit
           indCredit          - indices of forward problems to work on
           dumpResults        - A function pointer for saving the results throughout the iterations.
                                 - We assume that dumpResults is dumpResults(mc,Dc,iter,pInv,pMis),
                                 - where mc is the recovered model, Dc is the predicted data.
                                 - If dumpResults is not given, nothing is done (dummy() is called).

    Output:
           mc                 - final model
           Dc                 - final computed data
           tikhonovFlag       - Data misfit convergence flag
           hist               - Iteration history. This is a vector. Each entry is
                                 - a structure containing the projGNCG history for each
                                 - alpha value.
"""
function iteratedTikhonov(mc,pInv::InverseParam,pMis,alphaParam,
                          targetMisfit;indCredit=[],dumpResults::Function=dummy,
                          solveGN=projPCG)

    alphaMax = alphaParam[1]; alphaMin = alphaParam[2]
    alphaFac = alphaParam[3]; nAlpha   = alphaParam[4]
    iter = 0
    tikhonovFlag = -1
    hist = GNhis[]
    Dc = []
    while iter < nAlpha
        iter += 1
        println("Starting projGNCG minimization with alpha $iter of $nAlpha")
        println("alpha = $(pInv.alpha)")
        mc,Dc,GNflag,GNhist = projGN(mc,pInv,pMis,dumpResults=dumpResults)
        push!(hist,GNhist)
        pInv.mref  = mc
        pInv.alpha = pInv.alpha/alphaFac
        GNiterIdx  = findall(hist[end].F .> 0.0)
        println(hist[end].F)
        if minimum(hist[end].F[GNiterIdx]) <= targetMisfit
            tikhonovFlag = 1
            break
        end
    end

    if tikhonovFlag < 0
        println("iteratedTikhonov iterated maximum number of times, using $nAlpha alpha values but reached only a misfit of $(hist[end].F[end]).")
    else
        println("iteratedTikhonov exiting after reaching desired misfit")
    end
    return mc,Dc,tikhonovFlag,hist
end


text:
# FT and NaN test
println("\nTesting the structures...")
@testset "Hydraulics --- structures" begin
    for FT in [Float32, Float64]
        leaf  = LeafHydraulics{FT}();
        root  = RootHydraulics{FT}();
        stem  = StemHydraulics{FT}();
        grass = create_grass(FT(-2.1), FT(0.5), FT[0,-1,-2,-3], collect(FT,0:1:20));
        palm  = create_palm(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        tree  = create_tree(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        treet = TreeSimple{FT}();
        _vc1  = WeibullSingle{FT}();
        _vc2  = WeibullDual{FT}();
        _sh1  = BrooksCorey{FT}();
        _sh2  = VanGenuchten{FT}();

        # Test the struct
        for data_set in [ leaf, root, stem, grass, palm, tree, treet, _vc1, _vc2, _sh1, _sh2]
            @test FT_test(data_set, FT)
            @test NaN_test(data_set)
        end

        # Test Soil types
        _nams = ["Sand", "Loamy Sand", "Sandy Loam", "Loam", "Sandy Clay Loam",
                 "Silt Loam", "Silt", "Clay Loam", "Silty Clay Loam",
                 "Sandy Clay", "Silty Clay", "Clay"];
        for _name in _nams
            create_soil_VC(_sh1, _name);
            create_soil_VC(_sh2, _name);
        end
        println("Expecting a warning here");
        create_soil_VC(_sh1, "Other");
        @test true;
    end
end

# FT and NaN test
println("\nTesting the soil VC functions...")
@testset "Hydraulics --- soil VC" begin
    for FT in [Float32, Float64]
        _sh1  = BrooksCorey{FT}();
        _sh2  = VanGenuchten{FT}();

        # test soil functions
        _rwc = FT(0.5);
        _p   = FT(-0.5);
        for result in [ soil_erwc(_sh2, FT(1)),
                        soil_rwc(_sh1, _p),
                        soil_rwc(_sh2, _p),
                        soil_k_ratio_rwc(_sh1, _rwc),
                        soil_k_ratio_rwc(_sh2, _rwc),
                        soil_k_ratio_swc(_sh1, _rwc*_sh1.Θs),
                        soil_k_ratio_swc(_sh2, _rwc*_sh2.Θs),
                        soil_k_ratio_p25(_sh1, _p),
                        soil_k_ratio_p25(_sh2, _p),
                        soil_p_25_erwc(_sh1, FT(1)),
                        soil_p_25_erwc(_sh2, FT(1)),
                        soil_p_25_rwc(_sh1, _rwc),
                        soil_p_25_rwc(_sh2, _rwc),
                        soil_p_25_swc(_sh1, _rwc*_sh1.Θs),
                        soil_p_25_swc(_sh2, _rwc*_sh2.Θs) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # test if the soil_p_25 and soil_rwc converge
        _p_giv = FT(-2);
        _r_m_1 = soil_erwc(_sh1, _p_giv);
        _r_m_2 = soil_erwc(_sh2, _p_giv);
        _p_m_1 = soil_p_25_erwc(_sh1, _r_m_1);
        _p_m_2 = soil_p_25_erwc(_sh2, _r_m_2);
        @test _p_giv ≈ _p_m_1 ≈ _p_m_2;

        _r_m_1 = soil_rwc(_sh1, _p_giv);
        _r_m_2 = soil_rwc(_sh2, _p_giv);
        _p_m_1 = soil_p_25_rwc(_sh1, _r_m_1);
        _p_m_2 = soil_p_25_rwc(_sh2, _r_m_2);
        @test _p_giv ≈ _p_m_1 ≈ _p_m_2;

        _r_m_1 = soil_swc(_sh1, _p_giv);
        _r_m_2 = soil_swc(_sh2, _p_giv);
        _p_m_1 = soil_p_25_swc(_sh1, _r_m_1);
        _p_m_2 = soil_p_25_swc(_sh2, _r_m_2);
        @test _p_giv ≈ _p_m_1 ≈ _p_m_2;
    end
end

# FT and NaN test
println("\nTesting the xylem VC functions...")
@testset "Hydraulics --- xylem VC" begin
    for FT in [Float32, Float64]
        _vc1 = WeibullSingle{FT}();
        _vc2 = WeibullDual{FT}();
        _sh1 = BrooksCorey{FT}();
        _sh2 = VanGenuchten{FT}();
        _p   = FT(-0.5);
        _v   = FT(0.9);

        # Test xylem_p_crit
        f_st = FT(1)
        for result in [ xylem_p_crit(_vc1, f_st),
                        xylem_p_crit(_vc2, f_st) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # test the xylem_k_ratio
        for result in [ xylem_k_ratio(_vc1, _p, _v),
                        xylem_k_ratio(_vc2, _p, _v),
                        xylem_k_ratio(_vc1, FT(1), _v),
                        xylem_k_ratio(_vc2, FT(1), _v) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
    end
end

# FT and NaN test
println("\nTesting the leaf functions...")
@testset "Hydraulics --- leaf" begin
    for FT in [Float32, Float64]
        leaf  = LeafHydraulics{FT}();
        _f    = FT(0.01);

        # Test the struct
        _lr = xylem_risk(leaf, _f);
        _ec = critical_flow(leaf);

        for result in [_lr, _ec]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
    end
end

# FT and NaN test
println("\nTesting the legacy functions...")
@testset "Hydraulics --- legacy" begin
    for FT in [Float32, Float64]
        _ps   = FT(-0.1);
        _ft   = FT(1);
        _fsl  = FT(0.8);
        _fsh  = FT(0.2);
        _rsl  = FT(0.5);

        # Test the struct
        grass = create_grass(FT(-2.1), FT(0.5), FT[0,-1,-2,-3], collect(FT,0:1:20));
        palm  = create_palm(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        tree  = create_tree(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        treet = TreeSimple{FT}();
        pressure_profile!(treet, _ps, _ft);
        @test FT_test(treet, FT);
        @test NaN_test(treet);

        treet = TreeSimple{FT}();
        pressure_profile!(treet, _ps, _fsl, _fsh, _rsl);
        @test FT_test(treet, FT);
        @test NaN_test(treet);

        for _plant in [grass, palm, tree]
            pressure_profile!(_plant, SteadyStateMode(); update=true);
            @test FT_test(_plant, FT);
            @test NaN_test(_plant);
        end

        inititialize_legacy!(grass);
        inititialize_legacy!(palm);
        inititialize_legacy!(tree);
        inititialize_legacy!(treet);
        @test FT_test(treet, FT);
        @test NaN_test(treet);
    end
end

# FT and NaN test
println("\nTesting the temperature functions...")
@testset "Hydraulics --- temperature functions" begin
    for FT in [Float32, Float64]
        leaf  = LeafHydraulics{FT}();
        root  = RootHydraulics{FT}();
        treet = TreeSimple{FT}();
        grass = create_grass(FT(-2.1), FT(0.5), FT[0,-1,-2,-3], collect(FT,0:1:20));
        palm  = create_palm(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        tree  = create_tree(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        T     = rand(FT) + 298;

        # test the temperature functions
        temperature_effects!(treet);
        temperature_effects!(grass);
        temperature_effects!(palm );
        temperature_effects!(tree );
        temperature_effects!(leaf, T);
        temperature_effects!(root, T);

        for dataset in [leaf, root, treet, grass, palm, tree]
            @test FT_test(dataset, FT);
            @test NaN_test(dataset);
        end
    end
end

# FT and NaN test
println("\nTesting the root-related functions...")
@testset "Hydraulics --- root-related functions" begin
    for FT in [Float32, Float64]
        root  = RootHydraulics{FT}();
        grass = create_grass(FT(-2.1), FT(0.5), collect(FT,0:-0.5:-3.0), collect(FT,0:1:20));

        # test the root q functions
        _p1 = FT(0)
        _p2 = FT(-0.5)
        _p3 = FT(-1.0)
        for result in [ xylem_flow(root, _p1),
                        xylem_flow(root, _p2),
                        xylem_flow(root, _p3) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # test the root qs from q functions
        _f0 = FT(0);
        _f1 = FT(0.1);
        _f2 = FT(0.5);
        _ks = zeros(FT, 5);
        _ps = zeros(FT, 5);
        _qs = zeros(FT, 5);
        for _f in [_f0, _f1, _f2]
            roots_flow!(grass.roots, _ks, _ps, _qs, _f);
            @test FT_test(grass, FT);
            @test NaN_test(grass);
        end

        # test the recalculate_roots_flow function
        _ks = zeros(FT, 5);
        _ps = zeros(FT, 5);
        _qs = zeros(FT, 5);
        roots_flow!(grass.roots, _ks, _ps, _qs, FT(0.5));
        roots_flow!(grass, FT(0.5));
        @test NaN_test(grass);
    end
end

# FT and NaN test
println("\nTesting the pressure functions...")
@testset "Hydraulics --- pressure functions" begin
    for FT in [Float32, Float64]
        leaf  = LeafHydraulics{FT}();
        root  = RootHydraulics{FT}();
        stem  = StemHydraulics{FT}();
        grass = create_grass(FT(-2.1), FT(0.5), FT[0,-1,-2,-3], collect(FT,0:1:20));
        palm  = create_palm(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        tree  = create_tree(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        treet = TreeSimple{FT}();

        # test the end_pressure function
        _f_1 = FT(0.001)
        _f_2 = FT(1)
        for result in [ end_pressure(leaf, _f_1),
                        end_pressure(leaf, _f_2),
                        end_pressure(root, _f_1),
                        end_pressure(root, _f_2),
                        end_pressure(stem, _f_1),
                        end_pressure(stem, _f_2) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # test the end_pressure function for plants
        for result in [ end_pressure(treet, _f_1),
                        end_pressure(treet, _f_2),
                        end_pressure(treet, _f_1, _f_1, FT(0.5)) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
    end
end

# FT and NaN test
println("\nTesting the plant-level functions...")
@testset "Hydraulics --- plant-level" begin
    for FT in [Float32, Float64]
        grass = create_grass(FT(-2.1), FT(0.5), FT[0,-1,-2,-3], collect(FT,0:1:20));
        palm  = create_palm(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        tree  = create_tree(FT(-2.1), FT(5.5), FT(6), FT[0,-1,-2,-3], collect(FT,0:1:20));
        treet = TreeSimple{FT}();

        # test the critical_flow function
        for result in [ critical_flow(treet) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # test the plant conductances
        plant_conductances!(treet);
    end
end

# FT and NaN test
println("\nTesting the capacitance functions...")
@testset "Hydraulics --- capacitance" begin
    for FT in [Float32, Float64]
        grass = create_grass(FT(-2.1), FT(0.5), FT[0,-1,-2,-3], collect(FT,0:1:20));
        palm  = create_palm(FT(-2.1), FT(5.4), FT(8), FT[0,-1,-2,-3], collect(FT,0:1:20));
        tree  = create_tree(FT(-2.1), FT(5.4), FT(8), FT[0,-1,-2,-3], collect(FT,0:1:20));

        # test the critical_flow function
        update_PVF!(grass, FT(1)); @test true;
        update_PVF!(palm , FT(1)); @test true;
        update_PVF!(tree , FT(1)); @test true;
        for leaf in tree.leaves
            leaf.q_out = 2e-3;
        end
        update_PVF!(tree , FT(1)); @test true;
        update_PVF!(tree.roots[1], FT(1e6));
        update_PVF!(tree.roots[1], FT(1e6), true);
        update_PVF!(tree.trunk, FT(1e6));
        update_PVF!(tree.leaves[1], FT(1e6));

        # extra tests
        p_from_volume(PVCurveSegmented{FT}(), FT(0.9), FT(300));
        p_from_volume(PVCurveSegmented{FT}(), FT(0.6), FT(300));
        p_from_volume(PVCurveSegmented{FT}(), FT(0.1), FT(300));
    end
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module DiscreteAdjoint

using OrdinaryDiffEq, LinearAlgebra
using UnPack, MuladdMacro, PreallocationTools
import ForwardDiff, ReverseDiff, Zygote

using OrdinaryDiffEq: @.., True, False
using OrdinaryDiffEq: Tsit5Cache, Tsit5ConstantCache

include("derivatives.jl")
include("adjoint.jl")

include("caches/low_order_rk_caches.jl")

include("residuals/low_order_rk_residual.jl")

export discrete_adjoint
export ForwardDiffVJP, ReverseDiffVJP, ZygoteVJP

end


text:
## Calling Church from Julia
## =========================
function parse_output(out)
  lines = split(out,"\n")
  compile_time = lines[1][18:end]
  compile_time = int(filter(x->!(x == 'm' || x == 's'),compile_time))
  run_time = lines[2][18:end]
  run_time = int(filter(x->!(x == 'm' || x == 's'),run_time))
  data = lines[3]
  jl_data = split(data[2:end-1]," ")
  p_true = length(filter(s->s=="#t",jl_data))/length(jl_data)
  @compat Dict("compile_time" => compile_time, "run_time" => run_time, "prob" => p_true)
end

function run_church(program, mhsteps = 10:30:3000)
  [parse_output(readall(`church --timed --program-args $n /home/zenna/sandbox/$program`))
   for n = mhsteps]
end

function run_church(program, n::Integer)
  parse_output(readall(`church --timed --program-args $n /home/zenna/sandbox/$program`))
end


text:
module statmech_tm_solver

# Write your package code here.
__precompile__(true)

using LinearAlgebra
using Zygote
using Optim
using Random
using KrylovKit
using TensorKit
using TensorOperations
using JLD2
using Tullio
using LoopVectorization

using TensorKitAD
using ChainRules
using ChainRulesCore
import ChainRulesCore: rrule, frule
import Base: +, -, *, iterate, length, getindex, similar
import TensorKit: leftorth, rightorth
import LinearAlgebra: normalize, dot, ldiv!

# imps.jl
export  act,
        transf_mat,
        transf_mat_T,
        ovlp,
        nonherm_cost_func,
        lambda_gamma,
        left_canonical,
        iTEBD_truncate,
        ln_fidelity,
        free_energy,
        mps_add,
        left_canonical_QR,
        right_canonical_QR,
        entanglement_spectrum,
        entanglement_entropy,
        expand,
        tangent_map_tn,
        tangent_map,
        calculate_ALR,
        biorth_right_canonical,
        biorth_left_canonical

# utils.jl
export  toarray,
        get_chi,
        get_d,
        arr_to_TensorMap,
        convert_to_tensormap,
        quicksave,
        quickload,
        logsumexp

# mpo_zoo.jl
export  mpo_triangular_AF_ising,
        mpo_triangular_AF_ising_alternative,
        mpo_triangular_AF_ising_adapter,
        mpo_kink_processor,
        mpo_square_ising

# bi_direction.jl
export  bimps,
        A_canonical,
        B_canonical

# cmpo.jl
export  cmpo,
        get_phy,
        get_vir

# cmps.jl
export  cmps,
        get_matrices,
        convert_to_cmps,
        right_canonical,
        K_mat,
        log_ovlp,
        convergence_measure,
        finite_env,
        kinetic,
        particle_density,
        pairing,
        point_interaction,
        optimize_conv_meas,
        compress,
        truncation_check,
        gauge_fixing_proj,
        tangent_map,
        preconditioner,
        precond_grad,
        lieb_liniger_h_tangent_map

# cmpo_zoo.jl
export  cmpo_ising,
        cmpo_xy,
        cmpo_xz,
        energy_quantum_ising,
        cmpo_xxz,
        energy_quantum_xxz,
        cmpo_ising_realtime,
        energy_lieb_liniger,
        energy_pairing_boson

# bi_direction_q.jl
export  qbimps

# TensorKitAD_supp.jl
export  convert_to_array,
        elem_mult

include("imps.jl")
include("utils.jl")
include("mpo_zoo.jl")
include("bi_direction.jl")
include("cmpo.jl")
include("cmps.jl")
include("cmpo_zoo.jl")
include("bi_direction_q.jl")
include("TensorKitAD_supp.jl")

end


text:
const test_cases = [
    ("Empty", UInt8[], "text/plain; charset=utf-8"),
    ("Binary", UInt8[1, 2, 3], "application/octet-stream"),

    ("HTML document #1", b"blah blah blah", "text/html; charset=utf-8"),
    ("HTML document #2", b"", "text/html; charset=utf-8"),
    ("HTML document #3 (leading whitespace)", b"   <!DOCTYPE HTML>...", "text/html; charset=utf-8"),
    ("HTML document #4 (leading CRLF)", b"\r\n...", "text/html; charset=utf-8"),

    ("Plain text", b"This is not HTML. It has ☃ though.", "text/plain; charset=utf-8"),

    ("XML", b"\n<?xml!", "text/xml; charset=utf-8"),

    # Image types.
    ("GIF 87a", b"GIF87a", "image/gif"),
    ("GIF 89a", b"GIF89a...", "image/gif"),

    # Audio types.
    ("MIDI audio", UInt8['M','T','h','d',0x00,0x00,0x00,0x06,0x00,0x01], "audio/midi"),
    ("MP3 audio/MPEG audio", UInt8['I','D','3',0x03,0x00,0x00,0x00,0x00,0x0f], "audio/mpeg"),
    ("WAV audio #1", UInt8['R','I','F','F','b',0xb8,0x00,0x00,'W','A','V','E','f','m','t',' ',0x12,0x00,0x00,0x00,0x06], "audio/wave"),
    ("WAV audio #2", UInt8['R','I','F','F',',',0x00,0x00,0x00,'W','A','V','E','f','m','t',' ',0x12,0x00,0x00,0x00,0x06], "audio/wave"),
    ("AIFF audio #1", UInt8['F','O','R','M',0x00,0x00,0x00,0x00,'A','I','F','F','C','O','M','M',0x00,0x00,0x00,0x12,0x00,0x01,0x00,0x00,0x57,0x55,0x00,0x10,0x40,0x0d,0xf3,0x34], "audio/aiff"),
    ("OGG audio", UInt8['O','g','g','S',0x00,0x02,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x7e,0x46,0x00,0x00,0x00,0x00,0x00,0x00,0x1f,0xf6,0xb4,0xfc,0x01,0x1e,0x01,0x76,0x6f,0x72], "application/ogg"),

    # Video types.
    ("MP4 video", UInt8[0x00,0x00,0x00,0x18,'f','t','y','p','m','p','4','2',0x00,0x00,0x00,0x00,'m','p','4','2','i','s','o','m','<',0x06,'t',0xbf,'m','d','a','t'], "video/mp4"),
    ("AVI video #1", UInt8['R','I','F','F',',','O','\n',0x00,'A','V','I',' ','L','I','S','T','À'], "video/avi"),
    ("AVI video #2", UInt8['R','I','F','F',',','\n',0x00,0x00,'A','V','I',' ','L','I','S','T','À'], "video/avi"),
]

@testset "HTTP.sniff" begin
    for case in test_cases
        println("TEST - sniff.jl: ", case[1])
        @test HTTP.sniff(case[2]) == case[3]
    end

    @test HTTP.sniff(IOBuffer(test_cases[1][2])) == HTTP.sniff(test_cases[1][2])
end

json_strings = [
"n",
"nu",
"nul",
"null",
"t",
"tr",
"tru",
"true",
"f",
"fa",
"fal",
"fals",
"false",
"\"sample string \\\" with escaped double quote\"",
"1",
"1.0",
"3.1428",
"100000",
"[1,2,3]",
"[ \"simple array\" , {\"key\": null } ,{ \"key2\"  : false}]",
"{\"simpleobject\": null}",
"""
{
    "glossary": {
        "title": "example glossary",
		"GlossDiv": {
            "title": "S",
			"GlossList": {
                "GlossEntry": {
                    "ID": "SGML",
					"SortAs": "SGML",
					"GlossTerm": "Standard Generalized Markup Language",
					"Acronym": "SGML",
					"Abbrev": "ISO 8879:1986",
					"GlossDef": {
                        "para": "A meta-markup language, used to create markup languages such as DocBook.",
						"GlossSeeAlso": ["GML", "XML"]
                    },
					"GlossSee": "markup"
                }
            }
        }
    }
}
""",
"""
{"menu": {
  "id": "file",
  "value": "File",
  "popup": {
    "menuitem": [
      {"value": "New", "onclick": "CreateNewDoc()"},
      {"value": "Open", "onclick": "OpenDoc()"},
      {"value": "Close", "onclick": "CloseDoc()"}
    ]
  }
}}
""",
"""
{"widget": {
    "debug": "on",
    "window": {
        "title": "Sample Konfabulator Widget",
        "name": "main_window",
        "width": 500,
        "height": 500
    },
    "image": {
        "src": "Images/Sun.png",
        "name": "sun1",
        "hOffset": 250,
        "vOffset": 250,
        "alignment": "center"
    },
    "text": {
        "data": "Click Here",
        "size": 36,
        "style": "bold",
        "name": "text1",
        "hOffset": 250,
        "vOffset": 100,
        "alignment": "center",
        "onMouseUp": "sun1.opacity = (sun1.opacity / 100) * 90;"
    }
}}
""",
"""
{"web-app": {
  "servlet": [
    {
      "servlet-name": "cofaxCDS",
      "servlet-class": "org.cofax.cds.CDSServlet",
      "init-param": {
        "configGlossary:installationAt": "Philadelphia, PA",
        "configGlossary:adminEmail": "ksm@pobox.com",
        "configGlossary:poweredBy": "Cofax",
        "configGlossary:poweredByIcon": "/images/cofax.gif",
        "configGlossary:staticPath": "/content/static",
        "templateProcessorClass": "org.cofax.WysiwygTemplate",
        "templateLoaderClass": "org.cofax.FilesTemplateLoader",
        "templatePath": "templates",
        "templateOverridePath": "",
        "defaultListTemplate": "listTemplate.htm",
        "defaultFileTemplate": "articleTemplate.htm",
        "useJSP": false,
        "jspListTemplate": "listTemplate.jsp",
        "jspFileTemplate": "articleTemplate.jsp",
        "cachePackageTagsTrack": 200,
        "cachePackageTagsStore": 200,
        "cachePackageTagsRefresh": 60,
        "cacheTemplatesTrack": 100,
        "cacheTemplatesStore": 50,
        "cacheTemplatesRefresh": 15,
        "cachePagesTrack": 200,
        "cachePagesStore": 100,
        "cachePagesRefresh": 10,
        "cachePagesDirtyRead": 10,
        "searchEngineListTemplate": "forSearchEnginesList.htm",
        "searchEngineFileTemplate": "forSearchEngines.htm",
        "searchEngineRobotsDb": "WEB-INF/robots.db",
        "useDataStore": true,
        "dataStoreClass": "org.cofax.SqlDataStore",
        "redirectionClass": "org.cofax.SqlRedirection",
        "dataStoreName": "cofax",
        "dataStoreDriver": "com.microsoft.jdbc.sqlserver.SQLServerDriver",
        "dataStoreUrl": "jdbc:microsoft:sqlserver://LOCALHOST:1433;DatabaseName=goon",
        "dataStoreUser": "sa",
        "dataStorePassword": "dataStoreTestQuery",
        "dataStoreTestQuery": "SET NOCOUNT ON;select test='test';",
        "dataStoreLogFile": "/usr/local/tomcat/logs/datastore.log",
        "dataStoreInitConns": 10,
        "dataStoreMaxConns": 100,
        "dataStoreConnUsageLimit": 100,
        "dataStoreLogLevel": "debug",
        "maxUrlLength": 500}},
    {
      "servlet-name": "cofaxEmail",
      "servlet-class": "org.cofax.cds.EmailServlet",
      "init-param": {
      "mailHost": "mail1",
      "mailHostOverride": "mail2"}},
    {
      "servlet-name": "cofaxAdmin",
      "servlet-class": "org.cofax.cds.AdminServlet"},

    {
      "servlet-name": "fileServlet",
      "servlet-class": "org.cofax.cds.FileServlet"},
    {
      "servlet-name": "cofaxTools",
      "servlet-class": "org.cofax.cms.CofaxToolsServlet",
      "init-param": {
        "templatePath": "toolstemplates/",
        "log": 1,
        "logLocation": "/usr/local/tomcat/logs/CofaxTools.log",
        "logMaxSize": "",
        "dataLog": 1,
        "dataLogLocation": "/usr/local/tomcat/logs/dataLog.log",
        "dataLogMaxSize": "",
        "removePageCache": "/content/admin/remove?cache=pages&id=",
        "removeTemplateCache": "/content/admin/remove?cache=templates&id=",
        "fileTransferFolder": "/usr/local/tomcat/webapps/content/fileTransferFolder",
        "lookInContext": 1,
        "adminGroupID": 4,
        "betaServer": true}}],
  "servlet-mapping": {
    "cofaxCDS": "/",
    "cofaxEmail": "/cofaxutil/aemail/*",
    "cofaxAdmin": "/admin/*",
    "fileServlet": "/static/*",
    "cofaxTools": "/tools/*"},

  "taglib": {
    "taglib-uri": "cofax.tld",
    "taglib-location": "/WEB-INF/tlds/cofax.tld"}}}
""",
"""
{"menu": {
    "header": "SVG Viewer",
    "items": [
        {"id": "Open"},
        {"id": "OpenNew", "label": "Open New"},
        null,
        {"id": "ZoomIn", "label": "Zoom In"},
        {"id": "ZoomOut", "label": "Zoom Out"},
        {"id": "OriginalView", "label": "Original View"},
        null,
        {"id": "Quality"},
        {"id": "Pause"},
        {"id": "Mute"},
        null,
        {"id": "Find", "label": "Find..."},
        {"id": "FindAgain", "label": "Find Again"},
        {"id": "Copy"},
        {"id": "CopyAgain", "label": "Copy Again"},
        {"id": "CopySVG", "label": "Copy SVG"},
        {"id": "ViewSVG", "label": "View SVG"},
        {"id": "ViewSource", "label": "View Source"},
        {"id": "SaveAs", "label": "Save As"},
        null,
        {"id": "Help"},
        {"id": "About", "label": "About Adobe CVG Viewer..."}
    ]
}}
"""
]

@testset "HTTP.isjson" begin
    @test !HTTP.isjson("")[1]
    for str in json_strings
        @test HTTP.isjson(Vector{UInt8}(str))[1]
    end
end


text:
# common functionality

struct CompilerJob
    # core invocation
    f::Core.Function
    tt::DataType
    cap::VersionNumber
    kernel::Bool

    # optional properties
    minthreads::Union{Nothing,CuDim}
    maxthreads::Union{Nothing,CuDim}
    blocks_per_sm::Union{Nothing,Integer}
    maxregs::Union{Nothing,Integer}

    CompilerJob(f, tt, cap, kernel;
                    minthreads=nothing, maxthreads=nothing,
                    blocks_per_sm=nothing, maxregs=nothing) =
        new(f, tt, cap, kernel, minthreads, maxthreads, blocks_per_sm, maxregs)
end

# global job reference
# FIXME: thread through `job` everywhere (deadlocks the Julia compiler when doing so with
#        the LLVM passes in CUDAnative)
current_job = nothing

function signature(job::CompilerJob)
    fn = typeof(job.f).name.mt.name
    args = join(job.tt.parameters, ", ")
    return "$fn($(join(job.tt.parameters, ", ")))"
end

struct KernelError <: Exception
    job::CompilerJob
    message::String
    help::Union{Nothing,String}
    bt::StackTraces.StackTrace

    KernelError(job::CompilerJob, message::String, help=nothing;
                bt=StackTraces.StackTrace()) =
        new(job, message, help, bt)
end

function Base.showerror(io::IO, err::KernelError)
    println(io, "GPU compilation of $(signature(err.job)) failed")
    println(io, "KernelError: $(err.message)")
    println(io)
    println(io, something(err.help, "Try inspecting the generated code with any of the @device_code_... macros."))
    Base.show_backtrace(io, err.bt)
end

struct InternalCompilerError <: Exception
    job::CompilerJob
    message::String
    meta::Dict
    InternalCompilerError(job, message; kwargs...) = new(job, message, kwargs)
end

function Base.showerror(io::IO, err::InternalCompilerError)
    println(io, """CUDAnative.jl encountered an unexpected internal compiler error.
                   Please file an issue attaching the following information, including the backtrace,
                   as well as a reproducible example (if possible).""")

    println(io, "\nInternalCompilerError: $(err.message)")

    println(io, "\nCompiler invocation:")
    for field in fieldnames(CompilerJob)
        println(io, " - $field = $(repr(getfield(err.job, field)))")
    end

    if !isempty(err.meta)
        println(io, "\nAdditional information:")
        for (key,val) in err.meta
            println(io, " - $key = $(repr(val))")
        end
    end

    println(io, "\nInstalled packages:")
    for (pkg,ver) in Pkg.installed()
        println(io, " - $pkg = $ver")
    end

    println(io)
    versioninfo(io)
end

macro compiler_assert(ex, job, kwargs...)
    msg = "$ex, at $(__source__.file):$(__source__.line)"
    return :($(esc(ex)) ? $(nothing)
                        : throw(InternalCompilerError($(esc(job)), $msg;
                                                      $(map(esc, kwargs)...)))
            )
end

# maintain our own "global unique" suffix for disambiguating kernels
globalUnique = 0


text:
using BenchmarkTools
using Bobby

const SUITE = BenchmarkGroup()

SUITE["bitboard"] = BenchmarkGroup()
SUITE["bitboard"]["buildBoard"] = @benchmarkable Bobby.buildBoard()
SUITE["bitboard"]["buildLookUpTables"] = @benchmarkable Bobby.buildLookUpTables()

SUITE["moves"] = BenchmarkGroup()
b = Bobby.buildBoard()
l = Bobby.buildLookUpTables()
for f in (Bobby.getKingValid, Bobby.getNightsValid, Bobby.getBishopsValid,
	Bobby.getQueenValid, Bobby.getPawnsValid)
    SUITE["moves"][string(f)[7:end]] = @benchmarkable $(f)(b, l)
end
SUITE["moves"]["getRooksValid"] = @benchmarkable Bobby.getRooksValid(b)

SUITE["check"] = BenchmarkGroup()
SUITE["check"]["checkCheck"] = @benchmarkable Bobby.checkCheck(b)


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

function test_optmodattrs(uf, model, attr, listattr)
    @test !MOI.supports(model, attr)
    @test MOI.supports(uf, attr)
    @test isempty(MOI.get(uf, listattr))
    MOI.set(uf, attr, 0)
    @test MOI.get(uf, attr) == 0
    @test MOI.get(uf, listattr) == [attr]
end
function test_varconattrs(
    uf,
    model,
    attr,
    listattr,
    I::Type{<:MOI.Index},
    addfun,
    x,
    y,
    z,
)
    @test !MOI.supports(model, attr, I)
    @test MOI.supports(uf, attr, I)
    @test isempty(MOI.get(uf, listattr))
    MOI.set(uf, attr, [x, y], [2, 0])
    @test MOI.get(uf, attr, z) === nothing
    @test !MOI.is_empty(uf)
    @test MOI.get(uf, listattr) == [attr]
    MOI.set(uf, attr, z, 5)
    @test MOI.get(uf, attr, y) == 0
    @test MOI.get(uf, attr, [z, x]) == [5, 2]
    @test MOI.get(uf, listattr) == [attr]

    u = addfun(uf)
    @test MOI.get(uf, attr, u) === nothing
    @test MOI.get(uf, listattr) == [attr]
    MOI.set(uf, attr, u, 8)
    @test MOI.get(uf, listattr) == [attr]

    w = addfun(uf)
    @test MOI.get(uf, listattr) == [attr]
    @test MOI.get(uf, attr, w) === nothing

    @test MOI.is_valid(uf, u)
    MOI.delete(uf, u)
    @test !MOI.is_valid(uf, u)
    @test_throws MOI.InvalidIndex{typeof(u)} MOI.delete(uf, u)
    @test MOI.get(uf, listattr) == [attr]

    MOI.set(uf, attr, [w, z], [9, 4])
    @test MOI.get(uf, listattr) == [attr]
    @test MOI.get(uf, attr, w) == 9
    @test MOI.get(uf, attr, x) == 2
    @test MOI.get(uf, attr, z) == 4
    @test MOI.get(uf, attr, y) == 0
end

struct UnknownOptimizerAttribute <: MOI.AbstractOptimizerAttribute end

# A few objective/constraint types are supported to test both the fallback and the
# delegation to the internal model
MOIU.@model(
    ModelForUniversalFallback,
    (),
    (MOI.LessThan,),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    ()
)
function MOI.supports(
    ::ModelForUniversalFallback{T},
    ::Type{MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}},
) where {T}
    return false
end
function MOI.supports_constraint(
    ::ModelForUniversalFallback{T},
    ::Type{MOI.SingleVariable},
    ::Type{
        <:Union{
            MOI.EqualTo{T},
            MOI.GreaterThan{T},
            MOI.Interval{T},
            MOI.Integer,
            MOI.ZeroOne,
        },
    },
) where {T}
    return false
end

# TODO: Restructure to avoid sharing state across testsets. It doesn't seem
# necessary to use the same uf object for all the tests.
model = ModelForUniversalFallback{Float64}()
uf = MOIU.UniversalFallback(model)
@test MOI.is_empty(uf)
@testset "Copy Test" begin
    MOIT.copytest(uf, MOIU.Model{Float64}())
    @test !MOI.is_empty(uf)
    MOI.empty!(uf)
    @test MOI.is_empty(uf)
end
@testset "Start Values Test" begin
    src = MOIU.UniversalFallback(MOIU.Model{Float64}())
    dest = MOIU.UniversalFallback(MOIU.Model{Float64}())
    MOIT.start_values_test(dest, src)
end
@testset "Valid Test" begin
    MOIT.validtest(uf)
    @test !MOI.is_empty(uf)
    MOI.empty!(uf)
    @test MOI.is_empty(uf)
end
@testset "Empty Test" begin
    MOIT.emptytest(uf)
    @test MOI.is_empty(uf)
end
@testset "Name Test" begin
    MOIT.nametest(uf)
    @test !MOI.is_empty(uf)
    MOI.empty!(uf)
    @test MOI.is_empty(uf)
end
@testset "Optimizer Attribute" begin
    attr = UnknownOptimizerAttribute()
    listattr = MOI.ListOfOptimizerAttributesSet()
    empty!(uf.optattr)
    test_optmodattrs(uf, model, attr, listattr)
end
@testset "Model Attribute" begin
    attr = MOIT.UnknownModelAttribute()
    listattr = MOI.ListOfModelAttributesSet()
    test_optmodattrs(uf, model, attr, listattr)
    @test !MOI.is_empty(uf)
    MOI.empty!(uf)
    @test MOI.is_empty(uf)
end
x = MOI.add_variable(uf)
y, z = MOI.add_variables(uf, 2)
@testset "Variable Attribute" begin
    VI = MOI.VariableIndex
    attr = MOIT.UnknownVariableAttribute()
    listattr = MOI.ListOfVariableAttributesSet()
    test_varconattrs(uf, model, attr, listattr, VI, MOI.add_variable, x, y, z)
end
@testset "Objective Attribute" begin
    global x, y, z
    _single(vi) = MOI.SingleVariable(vi)
    function _affine(vi)
        return convert(
            MOI.ScalarAffineFunction{Float64},
            MOI.SingleVariable(vi),
        )
    end
    function _add_single_objective(vi)
        return MOI.set(
            uf,
            MOI.ObjectiveFunction{MOI.SingleVariable}(),
            _single(vi),
        )
    end
    function _add_affine_objective(vi)
        return MOI.set(
            uf,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            _affine(vi),
        )
    end
    @testset "Supported objective" begin
        F = MOI.SingleVariable
        @test MOI.supports(model, MOI.ObjectiveFunction{F}())
        _add_single_objective(x)
        @test MOI.get(uf, MOI.ObjectiveFunction{F}()) ≈ _single(x)
        MOI.set(uf, MOI.ObjectiveFunction{F}(), _single(y))
        @test MOI.get(uf, MOI.ObjectiveFunction{F}()) ≈ _single(y)
    end
    @testset "Unsupported objective" begin
        F = MOI.ScalarAffineFunction{Float64}
        @test !MOI.supports(model, MOI.ObjectiveFunction{F})
        @test MOI.supports(uf, MOI.ObjectiveFunction{F}())
        _add_affine_objective(x)
        @test MOI.get(uf, MOI.ObjectiveFunction{F}()) ≈ _affine(x)
        MOI.modify(
            uf,
            MOI.ObjectiveFunction{F}(),
            MOI.ScalarCoefficientChange(y, 1.0),
        )
        fx = MOI.SingleVariable(x)
        fy = MOI.SingleVariable(y)
        obj = 1fx + 1fy
        @test MOI.get(uf, MOI.ObjectiveFunction{F}()) ≈ obj
    end
end
@testset "Constraint Attribute" begin
    global x, y, z
    function _affine(vi)
        return convert(
            MOI.ScalarAffineFunction{Float64},
            MOI.SingleVariable(vi),
        )
    end
    function _add_affine_constraint(vi, ub)
        return MOI.add_constraint(uf, _affine(vi), MOI.LessThan(ub))
    end
    attr = MOIT.UnknownConstraintAttribute()
    @testset "Supported constraint" begin
        cx = _add_affine_constraint(x, 0.0)
        cy = _add_affine_constraint(y, 1.0)
        cz = _add_affine_constraint(z, 2.0)
        CI = MOI.ConstraintIndex{MOI.SingleVariable,MOI.LessThan{Float64}}
        listattr = MOI.ListOfConstraintAttributesSet{
            MOI.SingleVariable,
            MOI.LessThan{Float64},
        }()
        test_varconattrs(
            uf,
            model,
            attr,
            listattr,
            CI,
            uf -> _add_affine_constraint(x, 0.0),
            cx,
            cy,
            cz,
        )

        MOI.set(uf, MOI.ConstraintFunction(), cx, _affine(y))
        @test MOI.get(uf, MOI.ConstraintFunction(), cx) ≈ _affine(y)
        @test MOI.get(uf, MOI.CanonicalConstraintFunction(), cx) ≈ _affine(y)
        @test MOIU.is_canonical(
            MOI.get(uf, MOI.CanonicalConstraintFunction(), cx),
        )

        @test MOI.supports(uf, MOI.ConstraintName(), typeof(cx))
        MOI.set(uf, MOI.ConstraintName(), cx, "LessThan")
        @test MOI.get(uf, MOI.ConstraintName(), cx) == "LessThan"
        @test MOI.get(uf, typeof(cx), "LessThan") == cx
        MOI.delete(uf, cx)
        @test MOI.get(uf, typeof(cx), "LessThan") === nothing
    end
    # To remove the constraint attributes added in the previous testset
    MOI.empty!(uf)
    x = MOI.add_variable(uf)
    y, z = MOI.add_variables(uf, 2)
    @testset "Unsupported constraint" begin
        cx = _add_affine_constraint(x, 0.0)
        cy = _add_affine_constraint(y, 1.0)
        cz = _add_affine_constraint(z, 2.0)
        CI = MOI.ConstraintIndex{MOI.SingleVariable,MOI.EqualTo{Float64}}
        listattr = MOI.ListOfConstraintAttributesSet{
            MOI.SingleVariable,
            MOI.EqualTo{Float64},
        }()
        test_varconattrs(
            uf,
            model,
            attr,
            listattr,
            CI,
            uf -> _add_affine_constraint(x, 0.0),
            cx,
            cy,
            cz,
        )

        MOI.set(uf, MOI.ConstraintFunction(), cx, _affine(y))
        @test MOI.get(uf, MOI.ConstraintFunction(), cx) ≈ _affine(y)
        @test MOI.get(uf, MOI.CanonicalConstraintFunction(), cx) ≈ _affine(y)
        @test MOIU.is_canonical(
            MOI.get(uf, MOI.CanonicalConstraintFunction(), cx),
        )

        @test MOI.supports(uf, MOI.ConstraintName(), typeof(cx))
        MOI.set(uf, MOI.ConstraintName(), cx, "EqualTo")
        @test MOI.get(uf, MOI.ConstraintName(), cx) == "EqualTo"
        @test MOI.get(uf, typeof(cx), "EqualTo") == cx
        MOI.delete(uf, cx)
        @test MOI.get(uf, typeof(cx), "EqualTo") === nothing
    end
end
config = MOIT.TestConfig(solve = false)
@testset "empty" begin
    MOI.empty!(uf)
    @test MOI.is_empty(uf)
end
@testset "Unit" begin
    MOIT.unittest(uf, config)
end
@testset "Modification" begin
    MOIT.modificationtest(uf, config)
end
@testset "Continuous Linear" begin
    MOIT.contlineartest(uf, config)
end

@testset "Duplicate names" begin
    model = ModelForUniversalFallback{Float64}()
    uf = MOIU.UniversalFallback(model)
    x = MOI.add_variable(uf)
    # These are both intercepted by the fallback.
    function _affine(vi)
        return convert(
            MOI.ScalarAffineFunction{Float64},
            MOI.SingleVariable(vi),
        )
    end
    x_eq = MOI.add_constraint(uf, _affine(x), MOI.EqualTo(1.0))
    x_gt = MOI.add_constraint(uf, _affine(x), MOI.GreaterThan(1.0))
    MOI.set(uf, MOI.ConstraintName(), x_eq, "a name")
    MOI.set(uf, MOI.ConstraintName(), x_gt, "a name")
    @test_throws Exception MOI.get(
        uf,
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        },
        "a name",
    )
    @test_throws Exception MOI.get(uf, MOI.ConstraintIndex, "a name")
end

@testset "Deterministic constraint ordering" begin
    F = MOI.ScalarAffineFunction{Float64}
    _affine(vi) = convert(F, MOI.SingleVariable(vi))
    set1 = MOI.EqualTo(1.0)
    set2 = MOI.GreaterThan(1.0)
    for sets in [[set1, set2], [set2, set1]]
        model = ModelForUniversalFallback{Float64}()
        uf = MOIU.UniversalFallback(model)
        x = MOI.add_variable(uf)
        y = MOI.add_variable(uf)
        cx1 = MOI.add_constraint(uf, _affine(x), sets[1])
        cx2 = MOI.add_constraint(uf, _affine(x), sets[2])
        cy1 = MOI.add_constraint(uf, _affine(y), sets[1])
        cy2 = MOI.add_constraint(uf, _affine(y), sets[2])
        # check that the constraint types are in the order they were added in
        @test MOI.get(uf, MOI.ListOfConstraints()) ==
              [(F, typeof(sets[1])), (F, typeof(sets[2]))]
        # check that the constraints given the constraint type are in the order they were added in
        for set in sets
            @test MOI.get(uf, MOI.ListOfConstraintIndices{F,typeof(set)}()) == [
                MOI.ConstraintIndex{F,typeof(set)}(1),
                MOI.ConstraintIndex{F,typeof(set)}(2),
            ]
        end
    end
end

@testset "Delete" begin
    model = ModelForUniversalFallback{Float64}()
    uf = MOIU.UniversalFallback(model)
    MOIT.delete_test(uf)
end

@testset "Show" begin
    model = ModelForUniversalFallback{Float64}()
    uf = MOIU.UniversalFallback(model)
    @test sprint(show, uf) == MOI.Utilities.replace_acronym("""
    $(MOIU.UniversalFallback{ModelForUniversalFallback{Float64}})
    fallback for $(ModelForUniversalFallback{Float64})""")
end


text:
abstract type AbstractFuture end

abstract type AbstractSample end
abstract type AbstractSampleWithKeys <: AbstractSample end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
module SKCrossValidators

using PyCall

# standard included modules
using DataFrames
using Random
using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils

import AutoMLPipeline.CrossValidators.crossvalidate
export crossvalidate

function __init__()
    global SKM = pyimport("sklearn.metrics")
    global metric_dict = Dict(
          "roc_auc_score" => SKM.roc_auc_score,
          "accuracy_score" => SKM.accuracy_score,
          "auc" => SKM.auc,
          "average_precision_score" => SKM.average_precision_score,
          "balanced_accuracy_score" => SKM.balanced_accuracy_score,
          "brier_score_loss" => SKM.brier_score_loss,
          "classification_report" => SKM.classification_report,
          "cohen_kappa_score" => SKM.cohen_kappa_score,
          "confusion_matrix" => SKM.confusion_matrix,
          "f1_score" => SKM.f1_score,
          "fbeta_score" => SKM.fbeta_score,
          "hamming_loss" => SKM.hamming_loss,
          "hinge_loss" => SKM.hinge_loss,
          "log_loss" => SKM.log_loss,
          "matthews_corrcoef" => SKM.matthews_corrcoef,
          "multilabel_confusion_matrix" => SKM.multilabel_confusion_matrix,
          "precision_recall_curve" => SKM.precision_recall_curve,
          "precision_recall_fscore_support" => SKM.precision_recall_fscore_support,
          "precision_score" => SKM.precision_score,
          "recall_score" => SKM.recall_score,
          "roc_auc_score" => SKM.roc_auc_score,
          "roc_curve" => SKM.roc_curve,
          "jaccard_score" => SKM.jaccard_score,
          "zero_one_loss" => SKM.zero_one_loss
         )
end

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
                       sfunc::String,nfolds=10)
    @assert sfunc in keys(metric_dict)
    pfunc = metric_dict[sfunc]
    metric(X,Y) = pfunc(X,Y)
    crossvalidate(pl,X,Y,metric,nfolds)
end

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
                       sfunc::String,averagetype::String,nfolds=10)
    @assert sfunc in keys(metric_dict)
    pfunc = metric_dict[sfunc]
    metric(X,Y) = pfunc(X,Y,average=averagetype)
    crossvalidate(pl,X,Y,metric,nfolds)
end

end


text:
"""
```
compute_scenario_system(m, scen::Scenario; apply_altpolicy = false)
```

Given the current model parameters, compute the state-space system corresponding
to model `m` and alternative scenario `scen`. This function differs from
`compute_system` in that the `CCC`, `DD`, and `DD_pseudo` vectors are set to
zero (since we forecast in deviations from baseline) and shocks that are not in
`scen.instrument_names` are zeroed out in the `QQ` matrix.
"""
function compute_scenario_system(m::AbstractDSGEModel, scen::Scenario;
                                 apply_altpolicy::Bool = false)

    system = compute_system(m, apply_altpolicy = apply_altpolicy)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    # Zero out non-instrument shocks
    system.measurement.QQ = copy(system[:QQ])
    for shock in keys(m.exogenous_shocks)
        if !(shock in scen.instrument_names)
            shock_index = m.exogenous_shocks[shock]
            system[:QQ][shock_index, :] .= 0
            system[:QQ][:, shock_index] .= 0
        end
    end

    # Error out if there is nonzero measurement error
    if any(x -> x != 0, system[:EE])
        error("Can't simulate scenarios under nonzero measurement error")
    end

    return system
end

"""
```
filter_shocks!(m, scen, system::Scenario)
```

Given a scenario draw `scen`, back out the shocks necessary to hit
`scen.targets` and put them into `scen.instruments`. This function returns
`forecastshocks`, an `nshocks` x `horizon` matrix of filtered and smoothed
shocks.

This function checks `forecast_uncertainty_override(m)` for whether to smooth
shocks using the simulation smoother.
"""
function filter_shocks!(m::AbstractDSGEModel, scen::Scenario, system::System, in_sample::Bool = false)
    # Check applicability of this methodology
    @assert n_instruments(scen) >= n_targets(scen) "Number of instruments must be at least number of targets"

    # Construct data
    df = targets_to_data(m, scen)

    # Set initial state and state covariance to 0
    s_0 = zeros(n_states_augmented(m))
    P_0 = zeros(n_states_augmented(m), n_states_augmented(m))

    # Filter and smooth *deviations from baseline*
    _, forecastshocks, _ = smooth(m, df, system, s_0, P_0, draw_states = scen.draw_states,
                                  include_presample = true, in_sample = in_sample)

    # Assign shocks to instruments DataFrame
    for shock in keys(m.exogenous_shocks)
        shock_index = m.exogenous_shocks[shock]
        if shock in scen.instrument_names
            scen.instruments[shock] = forecastshocks[shock_index, :]
        else
            @assert all(x -> x == 0, forecastshocks[shock_index, :])
        end
    end

    return forecastshocks
end

"""
```
forecast_scenario_draw(m, scen::Scenario, system, draw_index)
```

Filter shocks and use them to forecast the `draw_index`th draw of `scen`.
"""
function forecast_scenario_draw(m::AbstractDSGEModel, scen::Scenario, system::System,
                                draw_index::Int)
    # Re-initialize model indices in case extra states or equations were added
    # for an alternative policy
    init_model_indices!(m)

    # Load targets
    load_scenario_targets!(m, scen, draw_index)

    # Filter shocks
    forecastshocks = filter_shocks!(m, scen, system)

    # Scale shocks if desired
    forecastshocks = scen.shock_scaling * forecastshocks

    # Re-solve model with alternative policy rule, if applicable
    if alternative_policy(m).key != :historical
        system = compute_scenario_system(m, scen, apply_altpolicy = true)
    end

    # Forecast
    s_T = zeros(n_states_augmented(m))
    forecaststates, forecastobs, forecastpseudo, _ =
        forecast(m, system, s_T, shocks = forecastshocks)

    # Check forecasted output matches targets *if not forecasting under
    # alternative policy or using simulation smoother*
    if alternative_policy(m).key == :historical && !scen.draw_states
        for var in scen.target_names
            var_index = m.observables[var]
            horizon = min(forecast_horizons(m), n_target_horizons(scen))
            @assert forecastobs[var_index, 1:horizon] ≈ scen.shock_scaling * scen.targets[1:horizon, var]
        end
    end

    # Return output dictionary
    forecast_output = Dict{Symbol, Array{Float64}}()
    forecast_output[:forecaststates] = forecaststates
    forecast_output[:forecastobs]    = forecastobs
    forecast_output[:forecastpseudo] = forecastpseudo
    forecast_output[:forecastshocks] = forecastshocks

    return forecast_output
end

"""
```
write_scenario_forecasts(m, scenario_output_files, forecast_output;
    verbose = :low)
```

Write scenario outputs in `forecast_output` to `values(scenario_output_files)`.
"""
function write_scenario_forecasts(m::AbstractDSGEModel,
                                  scenario_output_files::Dict{Symbol, String},
                                  forecast_output::Dict{Symbol, Array{Float64}};
                                  verbose::Symbol = :low)
    for (i, var) in enumerate([:forecastobs, :forecastpseudo])
        filepath = scenario_output_files[var]
        JLD2.jldopen(filepath, "w") do file
            write_forecast_metadata(m, file, var)
            write(file, "arr", Array{Float64}(forecast_output[var]))
            if :proportion_switched in keys(scenario_output_files)
                write(file, "proportion_switched", forecast_output[:proportion_switched][i])
            end
        end

        println(verbose, :high, " * Wrote " * basename(filepath))
    end
end

"""
```
forecast_scenario(m, scen::Scenario; verbose = :low)
```

Simulate all draws of `scen` using the modal parameters of the model `m`. This
function returns a `Dict{Symbol, Array{Float64}`.
"""
function forecast_scenario(m::AbstractDSGEModel, scen::Scenario;
                           verbose::Symbol = :low)
    # Print
    info_print(verbose, :low, "Forecasting scenario = " * string(scen.key) * "...")
    println(verbose, :low, "Start time: " * string(now()))
    println(verbose, :low, "Forecast outputs will be saved in " * rawpath(m, "scenarios"))

    start_time = time_ns()
    # Update model alt policy setting
    m <= Setting(:alternative_policy, scen.altpolicy, false, "apol",
                 "Alternative policy")

    # If no instrument names provided, use all shocks
    if isempty(scen.instrument_names)
        scen.instrument_names = collect(keys(m.exogenous_shocks))
    end

    # Load modal parameters and compute system
    params = load_draws(m, :mode; verbose = verbose)
    update!(m, params)
    system = compute_scenario_system(m, scen)

    # Get to work!
    ndraws = scen.n_draws == 0 ? count_scenario_draws!(m, scen) : scen.n_draws
    mapfcn = use_parallel_workers(m) ? pmap : map
    forecast_outputs = mapfcn(draw_ind -> forecast_scenario_draw(m, scen, system, draw_ind),
                              1:ndraws)

    # Assemble outputs and write to file
    forecast_outputs = convert(Vector{Dict{Symbol, Array{Float64}}}, forecast_outputs)
    forecast_output = assemble_block_outputs(forecast_outputs)
    output_files = get_scenario_output_files(m, scen, [:forecastobs, :forecastpseudo])
    write_scenario_forecasts(m, output_files, forecast_output, verbose = verbose)
    # Print
    forecast_time = (time_ns() - start_time)/1e9
    forecast_time_min = forecast_time/60
    println(verbose, :low, "\nTime elapsed: " * string(forecast_time_min) * " minutes")
    println(verbose, :low, "Forecast complete: " * string(now()))

    return forecast_output
end


text:
"All viable commands that can be sent to the engine"
Commands = (NORTH='n', SOUTH='s', EAST='e', WEST='w', STAY_STILL='o', GENERATE='g', CONSTRUCT='c', MOVE='m')

max_stars_repo_name: johnnykwwang/Halite-III

max_issues_repo_name: johnnykwwang/Halite-III

max_forks_repo_name: johnnykwwang/Halite-III


text:
module SparseQuadratureGrids

using JLD2
using Base.Cartesian
using StaticArrays

import Base:    length,
                getindex,
                setindex!,
                Val,
                convert

export  SplitWeights,
        QuadratureRule,
        KronrodPatterson,
        GenzKeister,
        NestedGrid,
        calc_Δ_prod!,
        smolyak!,
        FlatGrid,
        FlattenedGrid,
        GridContainer,
        GridBuild,
        AdaptiveBuild,
        aPrioriBuild,
        RawBuild,
        CacheBuild,
        Adaptive,
        AdaptiveRaw,
        Smolyak,
        SmolyakRaw,
        eval_grid!,
        calc_grid!,
        GridVessel,
        DynamicGridVessel,
        StaticGridVessel,
        default,
        SlidingVecFun,
        SlidingVec,
        SlidingVector,
        MatrixVecSVec,
        set!,
        eval!

#BLAS.set_num_threads(1)
include("univariate_rules.jl")
include("sliding_vectors.jl")
include("grid.jl")
include("grid_builders.jl")
include("grid_containers.jl")
include("grid_interfacing.jl")

end # module