Hello! I have a combinatorial discrete choice problem that I solve by Backward Induction. The problem seems to be intensive in both memory and computation time and I would be super happy to improve it on any of these two dimensions. Right now, the problem runs Out of Memory when K>12. I post a MWE
using BenchmarkTools
using ProgressMeter
using BitOperations
using Distributions
J = 1000
T = 20
Π = 100*rand(J,T)
Dᵢ = rand([0],J,1)
feat1 = rand([1,2,3,4,5,6,7,8],J,1)
feat2 = rand([1,2,3,4,5,6,7,8,9,10],J,1)
Φ1 = 10
Φ2 = 10
Sj = 20*rand(J,T)
Fj = 20*rand(J,T)
Π=Π.-Sj.-Fj
K = 10
Ind =([1,4,7,9,10,12,J-10,J-6,J-5,J-2])
# This is the Backward Induction Algorithm with interdependent choices
function myfun1(Ind,K,D::Array{Int64},Φ1,Φ2,feat1,feat2)
X1c = zeros(Int8,size(D,1))
X2c = zeros(Int8,size(D,1))
X1l = zeros(Int8,size(D,1))
X2l = zeros(Int8,size(D,1))
countc = zeros(Float64,size(D,1))
countl = zeros(Float64,size(D,1))
feat1=feat1[Ind]
feat2=feat2[Ind]
for j=1:K
for l=1:K
X1c[j] = X1c[j] + (feat1[j]==feat1[l])*D[l]
X1l[j] = X1l[j] + (feat2[j]==feat2[l])*D[l]
countc[j] = countc[j] + (feat1[j]==feat1[l])
countl[j] = countl[j] + (feat2[j]==feat2[l])
end
end
X1c = Φ1.*(X1c.>0)./countc
X1l = Φ2.*(X1l.>0)./countl
fcc = sum(X1c)
fcl = sum(X1l)
return fcc, fcl
end
function pshockprob(sequences,prob)
(KK,K)=size(sequences)
proba = zeros(KK)
for i=1:KK
proba[i]=prob^(sum(sequences[i,:]))*(1-prob)^(K-sum(sequences[i,:]))
end
return proba
end
function shocks2ind(sequences,shockpath)
(K,T)=size(shockpath)
inds = zeros(T)
shockpath[shockpath.==2] .= 0
for i=1:T
inds[i]=matchrow(shockpath[:,i]',sequences)
end
return convert(Array{Int64}, inds)
end
function myfun7(Ind,Π,Φ1,Φ2,Sj,prob,SeqShocks,Dᵢ,feat1,feat2,K,J,T)
β = 0.9
Π = Π[Ind,:]
Sj = Sj[Ind,:]
SeqShocks = SeqShocks[Ind,:]
D = zeros(Int64,2^K,2^K,K,T)
Indices = ones(Int64,2^K,2^K,T)
Dpast = ones(Int64,T+1)
v = zeros(Float64,2^K,2^K,T)
Dstar = zeros(Int64,K,T+1)
Dstar[:,1]=Dᵢ[Ind]
Πstar=0
ChoiceSet = zeros(Int64,2^K,K)
ShockSet = zeros(Int64,2^K,K)
GroupCosts = zeros(Float64,2^K,2^K,2)
let int_s = fill(false,K) #pre-allocate output
for j=0:2^K-1
@. int_s = bget(j,0:K-1)
ChoiceSet[j+1,:]=int_s
ShockSet[j+1,:]=int_s
end
end
for j=1:2^K, jj=1:2^K
(GroupCosts[j,jj,1],GroupCosts[j,jj,2])=myfun1(Ind,K,ChoiceSet[j,:].*ShockSet[jj,:],Φ1,Φ2,feat1,feat2)
end
pstatevec=pshockprob(ShockSet,prob)
indshocks=shocks2ind(ShockSet,SeqShocks)
#println(indshocks)
#println(ShockSet)
#println(pstatevec)
#println(GroupCosts)
vc = zeros(Float64, 2^K)
for t=0:T-1
for i=1:2^K, iii=1:2^K
#Basically what is wrong here is that I am still allowing people to enter
fill!(vc, 0.0)
for ii=1:2^K
for j=1:K
vc[ii] += ChoiceSet[ii,j]*(Π[j,T-t]+Sj[j,T-t]*ChoiceSet[i,j]-(1-ShockSet[iii,j])*999999999999999999999999999)
end
vc[ii] -= GroupCosts[ii,iii,1] + GroupCosts[ii,iii,2]
if t>0
vc[ii] += β*v[ii,:,T-t+1]'*pstatevec
end
end
#println(vc)
tempind = argmax(vc)
Indices[i,iii,T-t] = tempind
v[i,iii,T-t]=vc[tempind]
D[i,iii,:,T-t] = ChoiceSet[tempind,:]
end
end
for t=2:T+1
Dpast[t] = Indices[Dpast[t-1],indshocks[t-1],t-1]
Dstar[:,t]=D[Dpast[t],indshocks[t-1],:,t-1]
end
Dstar = Dstar[:,2:end]
return Dstar
end
println("TIMES")
p = 0.3
SeqShocks = (rand(J,T+1).>p).+1
BIC = @btime myfun7(Ind,Π,Φ1,Φ2,Sj,p,SeqShocks,Dᵢ,feat1,feat2,K,J,T)
println(sum(BIC))
Thanks so much in advance!