How can I transform this matlab code into Julia

Dear @skleinbo ,

how can I transform this matlab code into Julia code.

Here is the MATLAB code.

>> xx = zeros(2, 3)

xx =

     0     0     0
     0     0     0

>> yy = zeros(2, 3)

yy =

     0     0     0
     0     0     0

>> location = [1 2 3; 2 3 4]

location =

     1     2     3
     2     3     4

>> [xx(1,:), yy(1,:)] = ind2sub([2,2], location(1,:))

xx =

     1     2     1
     0     0     0


yy =

     1     1     2
     0     0     0

I want to transform this matlab code into Julia. and I try to rewrite the code,but get an error.

julia> xx = zeros(2, 3)
2×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> yy = zeros(2, 3)
2×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0


julia> location = [1 2 3; 2 3 4]
2×3 Matrix{Int64}:
 1  2  3
 2  3  4

julia> xx[1,:] .= CartesianIndices((2, 2))[location[1,:]][1]
ERROR: iteration is deliberately unsupported for CartesianIndex. Use `I` rather than `I...`, or use `Tuple(I)...`
Stacktrace:

Can you help me to correct my code to get the row and column index by the first row of location?

You can use LinearIndices and CartesianIndices to convert between linear and cartesian indices.

Edit: Ah, sorry, didn’t see you had discovered that already. You need to convert the cartesian indices to a tuple before indexing first, like the error suggests. Maybe

xx[1,:] .= first.(Tuple.(CartesianIndices((2, 2))[location[1,:]]))

As you can see, this is a bit unwieldy and inelegant and should prompt you to consider a less MATLABy way.

I try to use loop .

julia> for i in 1:2, j in 1:3
          xx[i, j] = CartesianIndices((2, 2))[location[i, j]][1]
       end

julia> for i in 1:2, j in 1:3
          yy[i, j] = CartesianIndices((2, 2))[location[i, j]][2]
       end

julia> xx
2×3 Matrix{Float64}:
 1.0  2.0  1.0
 2.0  1.0  2.0

julia> yy
2×3 Matrix{Float64}:
 1.0  1.0  2.0
 1.0  2.0  2.0

This is the result what I need. But I don’t think loop is the best way.

Thanks so much.This solves my problem.

Does Julia have the same function ind2sub or sub2ind? I think CartesianIndices or LinearIndices is worse than matlab function ind2sub.

Not for a long time PSA: replacement of ind2sub/sub2ind in Julia 0.7+

@skleinbo I have another problem. I have convert the results to Int, but why the type is still Float64?

 pol_ind_kp[iak_state, :] = Int.(first.(Tuple.(CartesianIndices((nk, na))[Int.(location[iak_state, :])])))
 
typeof(pol_ind_kp)
Matrix{Float64} (alias for Array{Float64, 2})

I have solved the problem by myself. I forget to set the zeros matrix Int when I define the pol_ind_kp.

Loops aren’t bad. They’re often the best way to solve a problem.

1 Like

Thanks for your reply. But I think loops may take a very long time.

They do in Matlab, but loops in Julia are fast.

2 Likes

Oh, It is really a good news to me. I use matlab and Stata very often. And loops always take a very long time. Thank you very much.

Matlab to Julia requires significant unlearning of bad habits. Another one for free: whenever you type rand(1) or rand(n,1), you almost certainly want rand() or rand(n).

Read carefully the Performance Tips · The Julia Language

Particularly the first one. There is where most commonly new Julia users fail to get a good performance.

Dear @skleinbo ,

how can I use LinearIndices to get the index. In matlab, we can use sub2ind.

>> A = [1 2 1; 2 2 1]

A =

     1     2     1
     2     2     1

>> inds = sub2ind([2,2], A(1,:), A(2,:))

inds =

     3     4     1

@ murphyk I try to use function sub2indv which @ murphyk wrote.

function subv2ind(shape, cindices)
    """Return linear indices given vector of cartesian indices.
    shape: d-tuple of dimensions.
    cindices: n-iterable of d-tuples, each containing cartesian indices.
    Returns: n-vector of linear indices.

    Based on:
    https://discourse.julialang.org/t/psa-replacement-of-ind2sub-sub2ind-in-julia-0-7/14666/8
    Similar to this matlab function:
    https://github.com/tminka/lightspeed/blob/master/subv2ind.m
    """
    lndx = LinearIndices(Dims(shape))
    n = length(cindices)
    out = Array{Int}(undef, n)
    for i = 1:n
        out[i] = lndx[cindices[i]...]
    end
    return out
end


julia> A = [1 2 1; 2 2 1]
2×3 Matrix{Int64}:
 1  2  1
 2  2  1

julia> subv2ind((2, 2), Tuple.(zip(A[1, :], A[2, :])))
3-element Vector{Int64}:
 3
 4
 1

Is there another simpler way to do this?

julia> subv2ind((2, 2), [(i, j) for (i, j) in zip(A[1, :], A[2, :])])
3-element Vector{Int64}:
 3
 4
 1
julia> sub2ind(shape, A) = map(eachcol(A)) do c
         LinearIndices(shape)[c...]
       end
sub2ind (generic function with 1 method)

julia> sub2ind((2,2), A)
3-element Vector{Int64}:
 3
 4
 1

This assumes that each column of A stores Cartesian indices. Works in higher dimensions too

julia> B = [1 2 1; 2 2 1; 2 1 2]
3×3 Matrix{Int64}:
 1  2  1
 2  2  1
 2  1  2

julia> sub2ind((2,2,2), B)
3-element Vector{Int64}:
 7
 4
 5

You seem to be translating MATLAB code one-to-one to Julia. This might not be the best idea. Julia is not MATLAB and doing things the “Julian” way is often more efficient, elegant, generic and faster. We can help you optimize a piece of Julia code.

@skleinbo Thanks a lot for your kind help. Can you help me to improve my Julia codes. It takes very long time. Here is the Julia codes.

using SpecialFunctions, LinearAlgebra #I需要使用LinearAlgebra,efrc需要使用SpecialFunctions
using XLSX,MATLAB
include("tauchen.jl")
include("ind_trans.jl")

x = 1

β = 0.87
δ = 0.06
γ = 2
α = 0.35
wage = 1
za = 1
ϵ = 3
μ = ϵ / (ϵ - 1)
D = 1

phil = 0.00

nzp = 2
zp_sd_target = 0.43
piL = 0.80
zLspace = range(0.0001, 0.9999, length=1000000)
difference = piL .* (zLspace .- 1) .^ 2 + (1 - piL) .* ((piL .- piL .* zLspace) ./ (1 - piL)) .^ 2 .- zp_sd_target^2;

value, index = findmin(abs.(difference))
zp_l = zLspace[index]

zp_h = (1 - piL * zp_l) / (1 - piL)
zp_grid = [zp_l, zp_h]

ntau = 1
if ntau == 2
    tau_grid = range(-0.29, 0.29, length=ntau)
    tauprob = [0.81 0.19; 0.19 0.81]
elseif ntau == 1
    tau_grid = [0]
    tauprob = 1
end

nq = 1
if nq == 2
    q_grid = range(-0.18, 0.18, length=nq)
elseif nq == 1
    q_grid = [0]
    qprob = 1
end


nzt = 11
rho_zt = 0.59
sigma_zt = 0.13
mu_zt = -(sigma_zt^2) / (2 * (1 + rho_zt))
@show zt_grid = tauchen(nzt, mu_zt, rho_zt, sigma_zt, 3)[1]
ztprob = tauchen(nzt, mu_zt, rho_zt, sigma_zt, 3)[2]

runexp = 1
nr = 6
r_grid = [0.01, 0.02, 0.03, 0.04, 0.05, 0.10]

if runexp == 0
    rho_r = 0.50
    sigma_r = 0.0086
    mu_r = 0.03 * rho_r
    r_grid = tauchen(nr - 1, mu_r, rho_r, sigma_r, 2.014)[1]
    rprob = tauchen(nr - 1, mu_r, rho_r, sigma_r, 2.014)[2]
    r_grid = [r_grid findmax(r_grid)[1]]
    rprob = [rprob zeros(nr - 1, 1); 0.00 0.00 0.00 0.00 0.00 1.00]
end

chi0 = 0.98
chi1 = 0.047
psi = 3.2
k_l = 0.01
k_h = 6.0
nk = 120
k_grid = range(k_l, k_h, length=nk)

a_l = 0.01
a_h = 3.0
na = 120
a_grid = range(a_l, a_h, length=na)

n_choice = nk * na
n_state = nzp * nzt * nr * ntau * nq

r_grid_ind = 1:nr
zt_grid_ind = 1:nzt
zp_grid_ind = 1:nzp
tau_grid_ind = 1:ntau
q_grid_ind = 1:nq

Q = repeat(q_grid, nzp * nzt * nr * ntau, 1)
Q_ind = repeat(q_grid_ind, nzp * nzt * nr * ntau, 1)


TAU = repeat(tau_grid, nq, 1)
TAU = sort!(TAU, dims=1) #dims=1按行排序 dims = 2按列排序

TAU = repeat(TAU, nzp * nzt * nr, 1)
TAU_ind = repeat(tau_grid_ind, nq, 1)
TAU_ind = sort!(TAU_ind, dims=1)
TAU_ind = repeat(TAU_ind, nzp * nzt * nr, 1)

R = repeat(r_grid, ntau * nq, 1)
R = sort!(R, dims=2)
R = repeat(R, nzp * nzt, 1)
R_ind = repeat(r_grid_ind, ntau * nq, 1)
R_ind = sort!(R_ind, dims=1)
R_ind = repeat(R_ind, nzp * nzt, 1)

ZT = repeat(zt_grid', nr * ntau * nq, 1)
ZT = sort!(ZT, dims=1)
ZT = repeat(ZT, nzp, 1)

ZT_ind = repeat(zt_grid_ind, nr * ntau * nq, 1)
ZT_ind = sort!(ZT_ind, dims=1)
ZT_ind = repeat(ZT_ind, nzp, 1)

ZP = repeat(zp_grid, nzt * nr * ntau * nq, 1)
ZP = sort!(ZP, dims=1)
ZP_ind = repeat(zp_grid_ind, nzt * nr * ntau * nq, 1)
ZP_ind = sort(ZP_ind, dims=1);

EXOG = [ZP ZT R TAU Q]
EXOG_ind = [ZP_ind ZT_ind R_ind TAU_ind Q_ind]
EXOG = EXOG'
EXOG_ind = EXOG_ind'

A = repeat(a_grid, nk, 1)
A = sort!(A, dims=1)
K = repeat(k_grid, na, 1)
G = [A K]

prob = zeros(n_state, n_state);

for i_state = 1:n_state
    izp = EXOG_ind[1, i_state]
    izt = EXOG_ind[2, i_state]
    ir = EXOG_ind[3, i_state]
    itau = EXOG_ind[4, i_state]
    iq = EXOG_ind[5, i_state]

    for i_state_next = 1:n_state
        izpnext = EXOG_ind[1, i_state_next]
        iztnext = EXOG_ind[2, i_state_next]
        irnext = EXOG_ind[3, i_state_next]
        itaunext = EXOG_ind[4, i_state_next]
        iqnext = EXOG_ind[5, i_state_next]

        prob[i_state, i_state_next] = ztprob[izt, iztnext] * I[izp, izpnext] * rprob[ir, irnext] * tauprob[itau, itaunext] * I[iq, iqnext]
    end
end


iter_vfi = 0
iter_howard = 0
error_vfi = 10^15
error_vfi_old = error_vfi
VFIcontinue = 1
error_tolerance = 10^(-6)

if γ < 2.5 && runexp == 1
    penalty = -10^(9)
else
    penalty = -10^(12)
end

Nh = 20

indicator = ones(n_choice, n_state)
pol_ind_kp = zeros(Int, n_choice, n_state)
pol_ind_ap = zeros(Int, n_choice, n_state)
V_new = zeros(n_choice, n_state)
V_old = penalty .* ones(n_choice, n_state)
Vh_old = penalty .* ones(n_choice, n_state)
Vh_new = penalty .* ones(n_choice, n_state);

collateral = chi0 .* G[:, 1] + chi1 .* (exp.(G[:, 2]) .- 1) .- G[:, 2]

for i_state = 1:n_state
    izp = EXOG_ind[1, i_state]
    izt = EXOG_ind[2, i_state]
    ir = EXOG_ind[3, i_state]
    itau = EXOG_ind[4, i_state]
    iq = EXOG_ind[5, i_state]

    for iak_state = 1:n_choice

        ik = Int(rem(iak_state - 1, nk) + 1)
        ia = Int(floor((iak_state - 1) / nk) + 1)

        dummy = ones(na * nk, 1)
        l = phil + ((za * exp(zt_grid[izt]) * zp_grid[izp])^((ϵ - 1) / (1 + α * (ϵ - 1)))) * (Complex(k_grid[ik] + q_grid[iq])^(α * (ϵ - 1) / (1 + α * (ϵ - 1)))) * (D^(ϵ / (1 + α * (ϵ - 1)))) * (μ^(-ϵ / (1 + α * (ϵ - 1)))) * ((wage * (1 + tau_grid[itau]) / (1 - α))^(-ϵ / (1 + α * (ϵ - 1))))
        y = (za * exp(zt_grid[izt]) * zp_grid[izp]) * (Complex(k_grid[ik] + q_grid[iq])^(α)) * ((l - phil)^(1 - α))
        T = wage * tau_grid[itau] * l + wage * (1 + tau_grid[itau]) * phil - (r_grid[ir] + δ) * q_grid[iq]
        pi = D * (y^((ϵ - 1) / ϵ)) - wage * (1 + tau_grid[itau]) * l

        c = pi - (r_grid[ir] + δ) * k_grid[ik] + (1 + r_grid[ir]) * a_grid[ia] + T .- G[:, 1] - (psi .* (G[:, 2] .- k_grid[ik]) .^ 2) ./ (2 * k_grid[ik])

        CAPeffective = k_grid[ik] + q_grid[iq]

        if CAPeffective <= 0
            dummy[:] .= -1
        end

        dummy[real.(c).<=0] .= -1
        dummy[collateral.<0] .= -1

        if findmax(dummy)[1] < 0
            indicator[iak_state, i_state] = 0
        end
    end
end

ZP_mat = repeat(transpose(EXOG[1, :]), na * nk, 1)
ZT_mat = repeat(transpose(EXOG[2, :]), na * nk, 1)
R_mat = repeat(transpose(EXOG[3, :]), na * nk, 1)
Q_mat = repeat(transpose(EXOG[5, :]), na * nk, 1)

AP_mat = repeat(G[:, 1], 1, n_state)
KP_mat = repeat(G[:, 2], 1, n_state)

l_vfi = zeros(Complex{Float64}, n_choice, n_state)
y_vfi = zeros(Complex{Float64}, n_choice, n_state)
T_vfi = zeros(Complex{Float64}, n_choice, n_state)
pi_vfi = zeros(Complex{Float64}, n_choice, n_state)

location = zeros(Int, n_choice, n_state)

EXP = zeros(n_choice, n_state);


for iak_state = 1:n_choice
    ik = CartesianIndices((nk, na))[iak_state][1]
    ia = CartesianIndices((nk, na))[iak_state][2]

    l_vfi[iak_state, :] = transpose(phil .+ ((za * exp.(EXOG[2, :]) .* EXOG[1, :]) .^ ((ϵ - 1) / (1 + α * (ϵ - 1)))) .* (Complex.(k_grid[ik] .+ EXOG[5, :]) .^ (α * (ϵ - 1) / (1 + α * (ϵ - 1)))) .* (D^(ϵ / (1 + α * (ϵ - 1)))) .* (μ^(-ϵ / (1 + α * (ϵ - 1)))) .* ((wage .* (1 .+ EXOG[4, :]) ./ (1 - α)) .^ (-ϵ / (1 + α * (ϵ - 1)))))
    y_vfi[iak_state, :] = transpose((za * exp.(EXOG[2, :]) .* EXOG[1, :]) .* (Complex.(k_grid[ik] .+ EXOG[5, :]) .^ (α)) .* ((l_vfi[iak_state, :] .- phil) .^ (1 - α)))
    T_vfi[iak_state, :] = transpose(wage .* EXOG[4, :] .* l_vfi[iak_state, :] .+ wage .* (1 .+ EXOG[4, :]) .* phil .- (EXOG[3, :] .+ δ) .* EXOG[5, :])
    pi_vfi[iak_state, :] = transpose(D .* (y_vfi[iak_state, :] .^ ((ϵ - 1) / ϵ)) - wage .* (1 .+ EXOG[4, :]) .* l_vfi[iak_state, :])
end

piplusT = zeros(na, n_choice)
c = zeros(na, n_choice)
V_sub = zeros(na, n_choice)
collateral = chi0 .* AP_mat + chi1 .* (exp.(KP_mat) .- 1) - KP_mat

while VFIcontinue == 1
    error_vfi_old = error_vfi
    for i_state in 1:n_state
        prob_mat = repeat(transpose(prob[i_state, :]), n_choice, 1)
        EXP[:, i_state] = sum(prob_mat[:, :] .* V_old[:, :], dims=2)
    end
    for iak_state in 1:n_choice
        ik = Int(rem(iak_state - 1, nk) + 1)
        ia = Int(floor((iak_state - 1) / nk) + 1)

        piplusT = repeat(transpose(pi_vfi[iak_state, :]), n_choice, 1) + repeat(transpose(T_vfi[iak_state, :]), n_choice, 1)
        c = piplusT .- (R_mat .+ δ) * k_grid[ik] + (1 .+ R_mat) * a_grid[ia] - AP_mat - (psi * (KP_mat .- k_grid[ik]) .* (KP_mat .- k_grid[ik])) ./ (2 * k_grid[ik])
        CAPeffective = KP_mat + Q_mat

        V_sub = ((1 ./ abs.(c) .^ ((γ - 1) * ones(n_choice, n_state))) .- 1) ./ (1 - γ) + β * EXP

        @. V_sub[(real(c)<=0)|(collateral<0)|(CAPeffective<=0)] = penalty
        V_new[iak_state, :] = findmax(V_sub, dims=1)[1]
        location[iak_state, :] = first.(Tuple.(findmax(V_sub, dims=1)[2]))

        pol_ind_kp[iak_state, :] = first.(Tuple.(CartesianIndices((nk, na))[location[iak_state, :]]))
        pol_ind_ap[iak_state, :] = last.(Tuple.(CartesianIndices((nk, na))[location[iak_state, :]]))

    end

    ### HOWARD's Improment Algorithm
    Vh_old = V_new
    iter_howard = 0

    while iter_howard <= Nh

        for iak_state = 1:n_choice
            ik = Int(rem(iak_state - 1, nk) + 1)
            ia = Int(floor((iak_state - 1) / nk) + 1)

            c = transpose(pi_vfi[iak_state, :] - (EXOG[3, :] .+ δ) * k_grid[ik] + (1 .+ EXOG[3, :]) * a_grid[ia] + T_vfi[iak_state, :] - a_grid[pol_ind_ap[iak_state, :]] - (psi * (k_grid[pol_ind_kp[iak_state, :]] .- k_grid[ik]) .* (k_grid[pol_ind_kp[iak_state, :]] .- k_grid[ik])) ./ (2 * k_grid[ik]))
            Vh_new[iak_state, :] = ((1 ./ abs.(c) .^ ((γ - 1) * ones(1, nzt * nzp * nr * ntau * nq))) .- 1) ./ (1 - γ) + β .* sum(prob .* Vh_old[location[iak_state, :], :], dims=2)
        end


        Vh_new = Vh_new .* indicator + Vh_old .* (1 .- indicator)
        Vh_old = copy(Vh_new)
        iter_howard += 1
    end

    ### HOWARD END ##


    if Nh == -1
        error_vfi = findmax(abs.(V_new[:] - V_old[:]))[1]
        V_old = copy(V_new)
    else
        error_vfi = findmax(abs.(Vh_new[:] - V_old[:]))[1]
        V_old = copy(Vh_new)
    end

    println("--------------------------------------")

    println("error_vfi_old = $error_vfi_old")
    println("error_vfi = $error_vfi")

    iter_vfi += 1

    println("iter_vfi = $iter_vfi")

    if error_vfi < error_tolerance
        VFIcontinue = 0
    end

end


V_final = zeros(na, nk, nzp, nzt, nr, ntau, nq)
a_prime_val = zeros(na, nk, nzp, nzt, nr, ntau, nq)
k_prime_val = zeros(na, nk, nzp, nzt, nr, ntau, nq)
c_val = zeros(Complex{Float64}, na, nk, nzp, nzt, nr, ntau, nq)
binding_collateral = zeros(na, nk, nzp, nzt, nr, ntau, nq)
sim_pol_ind_ap = zeros(na, nk, nzp, nzt, nr, ntau, nq)
sim_pol_ind_kp = zeros(na, nk, nzp, nzt, nr, ntau, nq)
effective_lambda = zeros(na, nk, nzp, nzt, nr, ntau, nq)

for i_state = 1:n_state
    izp = EXOG_ind[1, i_state]
    izt = EXOG_ind[2, i_state]
    ir = EXOG_ind[3, i_state]
    itau = EXOG_ind[4, i_state]
    iq = EXOG_ind[5, i_state]
    for iak_state = 1:n_choice
        ik = first(Tuple(CartesianIndices((nk, na))[iak_state]))
        ia = last(Tuple(CartesianIndices((nk, na))[iak_state]))

        V_final[ia, ik, izp, izt, ir, itau, iq] = V_new[iak_state, i_state]
        sim_pol_ind_ap[ia, ik, izp, izt, ir, itau, iq] = pol_ind_ap[iak_state, i_state]
        sim_pol_ind_kp[ia, ik, izp, izt, ir, itau, iq] = pol_ind_kp[iak_state, i_state]
        a_prime_val[ia, ik, izp, izt, ir, itau, iq] = a_grid[pol_ind_ap[iak_state, i_state]]
        k_prime_val[ia, ik, izp, izt, ir, itau, iq] = k_grid[pol_ind_kp[iak_state, i_state]]
        c_val[ia, ik, izp, izt, ir, itau, iq] = pi_vfi[iak_state, i_state] - (r_grid[ir] + δ) * k_grid[ik] + (1 + r_grid[ir]) * a_grid[ia] + T_vfi[iak_state, i_state] - a_prime_val[ia, ik, izp, izt, ir, itau, iq] - (psi * (k_prime_val[ia, ik, izp, izt, ir, itau, iq] - k_grid[ik]) .^ 2) ./ (2 * k_grid[ik])

        if pol_ind_kp[iak_state, i_state] < nk
            effective_lambda[ia, ik, izp, izt, ir, itau, iq] = chi0 + chi1 * (exp(k_grid[pol_ind_kp[iak_state, i_state]+1]) - 1) / a_prime_val[ia, ik, izp, izt, ir, itau, iq]
            if chi0 > 10^5 || chi1 > 10^5
                effective_lambda[ia, ik, izp, izt, ir, itau, iq] = 10^8
            end
            if chi0 * a_prime_val[ia, ik, izp, izt, ir, itau, iq] + chi1 * (exp(k_grid[pol_ind_kp[iak_state, i_state]+1]) - 1) - k_grid[pol_ind_kp[iak_state, i_state]+1] < 0
                binding_collateral[ia, ik, izp, izt, ir, itau, iq] = 1
            else
                binding_collateral[ia, ik, izp, izt, ir, itau, iq] = 0
            end
        elseif pol_ind_kp[iak_state, i_state] == nk
            binding_collateral[ia, ik, izp, izt, ir, itau, iq] = 0
        end

        if V_final[ia, ik, izp, izt, ir, itau, iq] == penalty
            V_final[ia, ik, izp, izt, ir, itau, iq] = NaN
            a_prime_val[ia, ik, izp, izt, ir, itau, iq] = NaN
            k_prime_val[ia, ik, izp, izt, ir, itau, iq] = NaN
            sim_pol_ind_ap[ia, ik, izp, izt, ir, itau, iq] = NaN
            sim_pol_ind_kp[ia, ik, izp, izt, ir, itau, iq] = NaN
            c_val[ia, ik, izp, izt, ir, itau, iq] = NaN
            binding_collateral[ia, ik, izp, izt, ir, itau, iq] = NaN
            effective_lambda[ia, ik, izp, izt, ir, itau, iq] = NaN
        end

    end
end

@skleinbo I am new to Julia.The codes above takes very long time, I am waiting for the results. Many thanks in advance.

Best,
Raymond