# General comments on my code for filling a matrix

**URL:** <https://discourse.julialang.org/t/general-comments-on-my-code-for-filling-a-matrix/100913>\
**Category:** General Usage\
**Created:** [June 27, 2023, 8:56pm UTC](https://discourse.julialang.org/t/general-comments-on-my-code-for-filling-a-matrix/100913 "2023-06-27T20:56:16Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Zerosum](https://avatars.discourse-cdn.com/v4/letter/z/b2d939/32.png) [@Zerosum](https://discourse.julialang.org/u/Zerosum)\
**Post date:** [June 27, 2023, 8:56pm UTC](https://discourse.julialang.org/t/general-comments-on-my-code-for-filling-a-matrix/100913/1 "2023-06-27T20:56:16Z")

</div>

I wondered if I could get comments on my code. I’ve reviewed the official performance tips and a few blogs discussing achieving high performance.

Unfortunately, since I have a very scant knowledge of computation and memory allocation, I couldn’t improve the following codes more. I understand that any advice will heavily depend on the specific context and structure, but I was wondering if I could get a brief/broad comment, for instance, on which part I could consider for faster speed.

The garbage collecting time is roughly 20%, and I did try profiling but honestly got lost on how to read the result. Any comments would be more than helpful and appreciated.

```julia
function tenure_prob(type::String, x::Real, dist_q::Vector{Float64}, surplus::Array{Float64,5})

    @unpack_Parameter param

    if type == "a"
        μ = μ_a
        σ = σ_a
    else 
        μ = μ_b
        σ = σ_b
    end
    
    # Preallocate distributions 
    pdf_ϵ = vcat(0.0,diff(cdf(Normal(0,σ_ϵ), grids.ϵ)))
    pdf_mq = vcat(0.0,diff(cdf(Normal(μ,σ), grids.mq)))

    # tenure 0 matrix
    output_0_success = zeros(grid_n, grid_m) 

    # tenure 1 matrix 
    output_1_success = zeros(grid_n, grid_m, grid_ϵ) 

    # tenure 2 matrix
    output_2_success = zeros(grid_n, grid_m, grid_ϵ, grid_ϵ)  

    # tenure 3 matrix
    output_3_success = zeros(grid_n, grid_m, grid_ϵ, grid_ϵ, grid_ϵ)  

    # tenure 4 matrix 
    output_4_success = zeros(grid_n, grid_m, grid_ϵ, grid_ϵ, grid_ϵ, grid_ϵ) 

    # tenure 5 matrix 
    output_5_success = zeros(grid_n, grid_m, grid_ϵ, grid_ϵ, grid_ϵ, grid_ϵ)  

    @inbounds @threads for (q,mq) in collect(product(1:grid_n,1:grid_m))

        ind_q = findfirst(surplus[:,1,x,1,mq] .> 0.0) # shouldn't matter y 
        ind_q = isnothing(findfirst(surplus[:,1,x,1,mq] .> 0.0)) ? grid_n+1 : ind_q 

        prob_0_success = (q >= ind_q) ? dist_q[q] * pdf_mq[mq] : 0
        output_0_success[q,mq] = prob_0_success

        ind_five_mq = findfirst(surplus[q,1,min(x+5,max_x+1),6,:] .> 0.0)
        ind_five_mq = isnothing(ind_five_mq) ? grid_m+1 : ind_five_mq

        prob_5_success = mq >= ind_five_mq

        val_mq = grids.mq[mq]

        ind_one_y = findfirst(surplus[q,:,min(x+1,max_x+1),2,mq] .> 0.0) 
        ind_two_y = findfirst(surplus[q,:,min(x+2,max_x+1),3,mq] .> 0.0)
        ind_three_y = findfirst(surplus[q,:,min(x+3,max_x+1),4,mq] .> 0.0)
        ind_four_y = findfirst(surplus[q,:,min(x+4,max_x+1),5,mq] .> 0.0)

        thr_one_y = isnothing(ind_one_y) ? 1e+8 : grids.y[ind_one_y] - val_mq

        for ϵ_1 in 1:grid_ϵ
            prob_1_success = grids.ϵ[ϵ_1] >= thr_one_y ? pdf_ϵ[ϵ_1] : 0 
            output_1_success[q,mq,ϵ_1] = prob_0_success * prob_1_success 

            thr_two_y = isnothing(ind_two_y) ? 1e+8 : 2 * (grids.y[ind_two_y] - val_mq) - grids.ϵ[ϵ_1]  

            for ϵ_2 in 1:grid_ϵ  
                prob_2_success = grids.ϵ[ϵ_2] >= thr_two_y ? pdf_ϵ[ϵ_2] : 0 
                output_2_success[q,mq,ϵ_1,ϵ_2] = output_1_success[q,mq,ϵ_1] * prob_2_success

                thr_three_y = isnothing(ind_three_y) ? 1e+8 : 3 * (grids.y[ind_three_y] - val_mq) - grids.ϵ[ϵ_1] - grids.ϵ[ϵ_2] 

                prob_3_success = ((grids.ϵ .>= thr_three_y) .* pdf_ϵ)
                @views output_3_success[q,mq,ϵ_1,ϵ_2,:] = output_2_success[q,mq,ϵ_1,ϵ_2] * prob_3_success
                        
                thr_four_y = [isnothing(ind_four_y) ? 1e+8 : 4 * (grids.y[ind_four_y] - val_mq) - grids.ϵ[ϵ_1] - grids.ϵ[ϵ_2] - val_ϵ for val_ϵ in grids.ϵ] 
                @views output_4_success[q,mq,ϵ_1,ϵ_2,:,:] = output_2_success[q,mq,ϵ_1,ϵ_2] * ((grids.ϵ .>= thr_three_y) .* pdf_ϵ) .* ((grids.ϵ' .>= thr_four_y) .* pdf_ϵ')                         
                @views output_5_success[q,mq,ϵ_1,ϵ_2,:,:] = output_4_success[q,mq,ϵ_1,ϵ_2,:,:] * prob_5_success 
            end
        end

    end

    return (prob_0_success = sum(output_0_success), 
    prob_1_success = sum(output_1_success), prob_1_fail = sum(output_0_success) - sum(output_1_success),
    prob_2_success = sum(output_2_success), prob_2_fail = sum(output_1_success) - sum(output_2_success),
    prob_3_success = sum(output_3_success), prob_3_fail = sum(output_2_success) - sum(output_3_success),
    prob_4_success = sum(output_4_success), prob_4_fail = sum(output_3_success) - sum(output_4_success),
    prob_5_success = sum(output_5_success), prob_5_fail = sum(output_4_success) - sum(output_5_success))
end

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 27, 2023, 9:20pm UTC](https://discourse.julialang.org/t/general-comments-on-my-code-for-filling-a-matrix/100913/2 "2023-06-27T21:20:26Z")

</div>

Your code is [not runnable](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757), it looks like it [uses global variables](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-untyped-global-variables) (`grids`), it allocates lots of small temporary arrays (e.g. `thr_four_y`, `prob_3_success`, `output_2_success[q,mq,ϵ_1,ϵ_2] * prob_3_success`, and many many others) in the inner loops, it assigns to slices (e.g. `output_4_success[q,mq,ϵ_1,ϵ_2,:,:]`) in the [wrong order for locality](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-column-major), and so on… there are so many problems it is hard to comment on all of them.

You’ll get more detailed advice if you make a _small_ (20 lines or so) _runnable_ example that illustrates a simplified version of what you are trying to do, instead of pasting in a large-ish function like this.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 27, 2023, 10:43pm UTC](https://discourse.julialang.org/t/general-comments-on-my-code-for-filling-a-matrix/100913/3 "2023-06-27T22:43:44Z")

</div>

> [@Zerosum](#):
>
> ```julia
> @inbounds @threads for (q,mq) in collect(product(1:grid_n,1:grid_m))
> 
> ```

This `@inbounds` will not apply to the `@threads` code, because `@threads` creates a closure, and `@inbounds` doesn’t penetrate it.

---

<div class="post-metadata">

**Author:** ![Zerosum](https://avatars.discourse-cdn.com/v4/letter/z/b2d939/32.png) [@Zerosum](https://discourse.julialang.org/u/Zerosum)\
**Post date:** [June 28, 2023, 6:15am UTC](https://discourse.julialang.org/t/general-comments-on-my-code-for-filling-a-matrix/100913/4 "2023-06-28T06:15:32Z")

</div>

Thank you very much!  
I’ve struggled to balance readability and giving up those small temporal arrays, but I see your point. I’ll try to develop an executable example that captures the core computation I’m doing.

---

<div class="post-metadata">

**Author:** ![Zerosum](https://avatars.discourse-cdn.com/v4/letter/z/b2d939/32.png) [@Zerosum](https://discourse.julialang.org/u/Zerosum)\
**Post date:** [June 28, 2023, 6:15am UTC](https://discourse.julialang.org/t/general-comments-on-my-code-for-filling-a-matrix/100913/5 "2023-06-28T06:15:58Z")

</div>

I didn’t know that. Thanks!