Collect values from function in Julia

Question

My question is regarding the use of the broadcast operator in Julia.

Suppose I have the following objects

M = [0.7 0.3; 0.4 0.6];
x0 = [100 100];
N=5;
y = zeros(N, size(x0)[2]);

function Markov_bling_bling(;Trans_mat, initial_states, n_ahead)
    # Define useful names
    M = Trans_mat; x0 = initial_states; N = n_ahead;
    # Compute the N-th state 
    xn = x0 * M^N
    return(x_n = xn)
end

(Sorry for the silly name)

So this function returns a 1x2 vector.

So I would like to store every xn as a row in y.

In R I would do this:

y <- list()
for(t in 1:(N+1)){
  y[t] = Markov_bling_bling(Trans_mat = M, initial_states = x0, n_ahead=(t-1))
}

y <- Reduce(rbind,x)

How can I accomplish this in Julia? I know that I have to use the broadcast operator in order to avoid a for loop.

But I still don't get how can I store the results, should I define y = []? What is Julia way to store results?

Thanks in advance!

Answer 1

It can be written in almost the same way as in R:

julia> reduce(vcat, [Markov_bling_bling(Trans_mat = M, initial_states = x0, n_ahead=(t-1)) for t in 1:N])
5×2 Matrix{Float64}:
 100.0   100.0
 110.0    90.0
 113.0    87.0
 113.9    86.1
 114.17   85.83

Answer 2

Bogumił's answer is the better direct translation, but may I suggest a more efficient version:

julia> function blingblings(x0, M, N)
           ys = similar(x0, size(x0, 2), N)
           ys[:, 1] = x0
           for i = 2:N
               ys[:, i] = ys[:, i-1]' * M
           end
           return ys
       end
blingblings (generic function with 1 method)

julia> blingblings(x0, M, 5)
2×5 Matrix{Float64}:
 100.0  110.0  113.0  113.9  114.17
 100.0   90.0   87.0   86.1   85.83

This, most importantly, avoids the quadratic complexity of exponentiating every time, and also respects the column-major nature of Julia arrays better. Nothing wrong with for-loops.

(I wanted to post a solution using accumulate, but unfortunately its init argument does not end up in the result, as I'd have expected, making that solution inelegant for your requirements:

julia> ys = accumulate(1:5; init=x0) do y, _
           (y * M)::Matrix{Float64}
       end
5-element Vector{Matrix{Float64}}:
 [110.0 90.0]
 [113.0 87.0]
 [113.89999999999999 86.1]
 [114.16999999999999 85.82999999999998]
 [114.25099999999998 85.74899999999998]

)