My goal is to calculate all the possible outcomes in a finite dynamic programming problem. The grid size is 100, there are three state variables and T periods.
A bottleneck is the following function, which I’ve started to improve in version 2. According to profiling, the issue is in line 28. Something about abstractarray.jl and multidimensional.jl. I’ve read that this is expected when indexing arrays. But I don’t know what can be done to speed up the code in spite of it. I’ve tried using inbounds and views, but these didn’t help much.
Also, according to code_warntype, Julia seems to have problems inferring the data types in lines 25 and 28.
Would you please have a look at it and tell me how I can improve it? Any help is appreciated.
Thanks!
using LinearAlgebra, Statistics
using BenchmarkTools
#packages to increase the speed of the code
using Profile # @profile command
using ProfileView #profiling using graphs
#number of periods
T = 13
#grid size
grid = 100
#continuation values
#period t (not t+1), θkt, ℓkt, ℓmt
s_W = zeros(T+1, grid, grid, grid)
j_W = similar(s_W)
#optimal payments
#period, θkt, ℓmt, outputs
Pst_array = zeros(T, grid, grid, 6) #we won't need payment proposal for period T, a.k.a. "T+1"
Pjt_array = similar(Pst_array)
#probability mass function
pmf = Array{Float64}(undef,grid,grid)
#example
pmf[:,:] .= 0.01
t, hkt, lkt, lmt, Pmt_array, k_W = T, 50, 50, 50, Pjt_array, s_W
continuation_value = 0.011
payment_threshold = 91
s_W[T+1, :, :, :] .= continuation_value
Pjt_array[t, payment_threshold:end, lkt, 1] .= continuation_value
#function to calculate the payment offer threshold ϕmt
#test with searchsorted first
function threshold_m(t, lkt, lmt, Pmt_array)
#searches the index where the payment offer is positive
ts = @views searchsortedfirst(Pmt_array[t, :, lkt, 1], 0.01)
#returns error message if threshold does not exist (101 is returned if missing by searchsortedfirst)
@assert ts < 101
return ts
end
#test
# t, lkt,lmt, Pmt_array = T, 50, 51, Pjt_array
# @btime threshold_m(t, lkt,lmt, Pmt_array)
# @code_warntype threshold_m(t, lkt,lmt, Pmt_array)
function respond_paymentv1(t, hkt, lkt, lmt, Pmt_array, k_W)
ts = threshold_m(t, lkt, lmt, Pmt_array)
lmt = max(ts, lmt)
#the possible values of payment given by the opponent
Pmt_grid = Pmt_array[t, lmt:end, lkt, 1]
#the screening cutoffs associated with each payment proposal, simplified
cutoffs = Array{Int64}(undef, length(Pmt_grid))
cutoffs .= lkt
#A1###
#payoff when accepting the payment is the best response
pA1 = copy(Pmt_grid)
#weighting by the probabilities
probk = pmf[hkt, hkt:end]
probm = pmf[lmt, lmt:end]
lowerbounds = collect(lmt:1:grid)
#K and M to make the shape of the arrays easier to see
K, M = length(probk), length(probm)
pA2 = Array{Float64}(undef, K, M)
for i in eachindex(cutoffs)
pA2[:, i] .= (k_W[t+1, hkt:end, cutoffs[i], lowerbounds[i]])
end
pA2
#indicator function
# IA1 = np.where(k_W[t+1, hkt:, cutoffs, np.arange(lmt, grid, step = 1)].T <= Pmt_grid, 1, 0) #python code
IA1 = Array{Int64}(undef, K, M)
IA1[findall(pA2 .<= transpose(repeat(Pmt_grid, 1, K)))] .= 1
#IA1 has S rows and J columns
#we want a matrix of probabilities with the same dimension
# prob = np.outer(probk, probm) #python code
# prob = repeat(probk, 1, M) .* transpose(repeat(probm, 1, K))
prob = probm' .* probk #outer product
pA1 = transpose(pA1) .* prob
A1 = pA1 .* IA1
#A2####
IA2 = @. ifelse(IA1==1, 0, 1)
pA2 = pA2 .* prob
A2 = pA2 .* IA2
payoff = A1 + A2
payoff = sum(payoff)
#returns the expected value of answering a payment offer and the probability of liquidation
return payoff
end
#test1: expected payoff of S
t, hkt, lkt, lmt, Pmt_array, k_W = T, 50, 50, 50, Pjt_array, s_W
@btime respond_paymentv1(t, hkt, lkt, lmt, Pmt_array, k_W)
function respond_paymentv2(t, hkt, lkt, lmt, Pmt_array, k_W)
ts = threshold_m(t, lkt, lmt, Pmt_array)
lmt = max(ts, lmt)
#the possible values of payment given by the opponent
Pmt_grid = Pmt_array[t, lmt:end, lkt, 1]
#the screening cutoffs associated with each payment proposal, simplified
cutoffs = Array{Int64}(undef, length(Pmt_grid))
cutoffs .= lkt
payoff = zeros(length(Pmt_grid))
# payoff = Array{Float64}(undef, length(Pmt_grid))
for (i, pmt) in enumerate(Pmt_grid)
#receives pmt if his skill level is below the cutoff
payoff[i] += pmt * sum(pmf[hkt, hkt:(cutoffs[i]-1)]) * pmf[lmt, lmt+(i-1)]
#receives pmt if his skill level greather than or equal to the cutoff
payoff[i] += dot(k_W[t+1, cutoffs[i]:end , cutoffs[i], lmt+(i-1)], pmf[hkt, cutoffs[i]:end]) * pmf[lmt, lmt+(i-1)]
end
return sum(payoff)
end
#test1: expected payoff of S
t, hkt, lkt, lmt, Pmt_array, k_W = T, 50, 50, 50, Pjt_array, s_W
@btime respond_paymentv2(t, hkt, lkt, lmt, Pmt_array, k_W)
@code_warntype respond_paymentv2(t, hkt, lkt, lmt, Pmt_array, k_W)
using Profile
Profile.clear()
for i in 1:10000
@profile respond_paymentv2(t, hkt, lkt, lmt, Pmt_array, k_W)
end
ProfileView.view()
PS1: Context
I’m studying a game where two players bargain with each other. I’ll call K the player who is taking the actions and M his adversary.
The gamel has asymmetric information. Each player bargaining power increases on his true skill level, hkt, which only he knows. But there are lower bounds, lkt and lmt, that help the adversary to infer hkt.
The function I’m struggling with calculates the expected payoff when K answers a payment proposal. K receives the payment pmt(hmt) in case he accepts it and k_W(hkt+1) if he rejects it. So he chooses the maximum between the two options. There are different possible values for hmt and for hkt+1, so we use the pmf (probability mass function) to calculate the expected payoff.
There is a cutoff level which is the minimum level of hkt+1 necessary to be lucrative to reject a payment offer. I’ve used this in version 2 to speed up the code.
PS2: If you have interest, the complete code is here. Again, tips are welcome