Hi!
I have a complementarity problem in which I am using the new awesome additions to JuMP to solve it. The problem involves a summation over 2^J terms, where J is around 30 and performance is critical.
I was wondering if I could get any recommendations on how to improve the performance of the code either on the proper JuMP function and on how to handle such a large problem.
Find a MWE below:
using JuMP, Combinatorics, BitOperations, BenchmarkTools
import PATHSolver
const J = 10
const param1 = 0.5
const param2 = 0.1
const prob = 0.5.*ones(J)
function possible_shock_combinations(J)
ShockSet = Array{Union{Nothing,Int}}(nothing, 2^J, J)
let int_s = fill(false, J) #pre-allocate output
for j = 0:2^J-1
@. int_s = bget(j, 0:J-1)
ShockSet[j+1, :] = int_s
end
end
return convert(Matrix{Int},ShockSet)
end
function joint_probability(p::Vector{Float64}, outcomes::Vector{Int})
# Check if the lengths of p and outcomes are the same
if length(p) != length(outcomes)
println("The lengths of probabilities and outcomes must be the same.")
end
prob = 1.0
for i in 1:length(outcomes)
if outcomes[i] == 1
prob *= p[i]
else
prob *= (1 - p[i])
end
end
return prob
end
function shockprobs(sequences,probs)
(KK, K) = size(sequences)
proba = zeros(Float64, KK)
for i = 1:KK
proba[i] = joint_probability(probs, sequences[i, :])
end
return proba
end
const shocks = possible_shock_combinations(J)
probs = shockprobs(shocks,prob) #probability of each event#
function Comp(param1, param2, shock, probs, w, i)
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, x[1:J] >= 0, start = 0.5)
@constraints(
model,
begin
c[j in 1:J], (-(1-param1)*w[i]*sum(probs[s]*(1-shock[s,j]*(1-param2))*sum(x[l]*(1 -shock[s,l]*(1-param2)) for l in 1: J)^(-1) for s in 1:2^J) + w[j])⟂ x[j]
end)
JuMP.optimize!(model)
return value.(x)
end
W = ones(J)
test = @btime Comp(param1, param2, shocks, probs, W, 1)