Which approach for multidimensional computation?

Specific Domains Data
,
1

Hello,
I need to make a multi-dimension computation, specifically

V_{ft,dc,r,t+1} = V_{ft,dc,r,t} *(1-m_{ft,dc,r,t} - 1/tp_{ft,dc,r,t}) + V_{ft_,dc-1,r,t}*tp_{ft,dc-1,r,t}*\beta_{ft,dc,r,t}

where, if you are curious, the contest is forest volumes that each year, for each forest time ft, region r and diameter class dc are equal to those present the previous year in the same cell, less those that die of mortality m and those that are “promoted” to the next diameter class (after time of passage tp) and plus those that instead promote from the previous diameter class last year (gaining \beta times their volume).

Now which is the best practice approach to develop this kind of model? I could use nested for loops - but that’s tedious, that in turn is one of the reason I would like to do it in Julia.
I could use matrix algebra (that I suppose would be the most efficient way), but I am not sure that my input data is all ordered the same… and the output too would not have labels associated to each dimension, that instead I would like to keep.
So, is there a way to use dataframes or dictionaries for this kind of multidimensional formulas ?
Which is the approach that you would suggest ?

2

Have you looked at packages which will generate the loops for you?

3

Would it work also for dictionaries? Because I am able to get it working only by position with arrays:

vol = Dict(("LO","hf",2015) => 1,
("LO","hf",2016) => 2,
("LO","con",2015) => 3,
("LO","con",2016) => 4,
("AQ","hf",2015) => 5,
("AQ","hf",2016) => 6,
("AQ","con",2015) => 7,
("AQ","con",2016) => 8,
)
m = Dict(("LO","hf",2015) => 0.11,
("LO","hf",2016) => 0.12,
("LO","con",2015) => 0.13,
("LO","con",2016) => 0.14,
("AQ","hf",2015) => 0.15,
("AQ","hf",2016) => 0.16,
("AQ","con",2015) => 0.17,
("AQ","con",2016) => 0.18,
)
@einsum deadV[r,ft,y] := vol[r,ft,y]* m[t,ft,y]
MethodError: no method matching size(::Dict{Tuple{String,String,Int64},Int64}, ::Int64)

vol = [a + b + c for a in 1.0:1.0:3.0, b in 3.0:1.0:4.0, c in 2.0:1.0:3.0]
m   = [a + b + c for a in 0.01:0.01:0.03, b in 0.02:0.01:0.03, c in 0.05:0.01:0.06]
@einsum deadV[r,ft,y] := vol[r,ft,y]* m[r,ft,y]
3×2×2 Array{Float64,3}:
[:, :, 1] =
 0.48  0.63
 0.63  0.8 
 0.8   0.99

[:, :, 2] =
 0.63  0.8 
 0.8   0.99
 0.99  1.2
4

Cartesian.jl is flexible: https://github.com/timholy/Cartesian.jl
EDIT:
Look at “Principles of Usage” for a brief overview. You can specify an anonymous function to specify the iterators that gives you a lot more flexibility than just 1:size(A,dim) used in most examples.

5

The modern version of that is Base.Cartesian

https://docs.julialang.org/en/latest/devdocs/cartesian/

6

Ah, oops. Using Base.Cartesian:

using Base.Cartesian
d1 = Dict("LO" => 4.2, "AQ" => 9.3)
d2 = Dict("hf" => 0.3, "con" => 111.1)
iter_strings_for_dim = (["LO", "AQ"],["hf","con"])
@nloops 2 i d -> begin
  iter_strings_for_dim[d]
end begin
  println(i_1, " ", i_2, ": ", d1[i_1] * d2[i_2])
end

Lets you iterate over strings.
Using the example dictionaries:

julia> iter_by_dim = (["LO", "AQ"],["hf","con"],[2015,2016])
(String["LO", "AQ"], String["hf", "con"], [2015, 2016])

julia> @nloops 3 i d -> begin
         iter_by_dim[d]
       end begin
         out += vol[i_1, i_2, i_3] * m[i_1, i_2, i_3]
       end

julia> out
5.640000000000001

There’s probably a better way to handle the different iteration per dimension.

7

uh… I am not sure that’s most concise/readable than just using list comprehension:

regions = ["LO","AQ"]
years = collect(2015:2016)
forTypes = ["hf","con"]
vol = Dict(("LO","hf",2015) => 100.0,
           ("LO","con",2015) => 300.0,
           ("AQ","hf",2015) => 500.0,
           ("AQ","con",2015) => 700.0,
)
m = Dict(("LO","hf",2015) => 0.11,
         ("LO","hf",2016) => 0.12,
         ("LO","con",2015) => 0.13,
         ("LO","con",2016) => 0.14,
         ("AQ","hf",2015) => 0.15,
         ("AQ","hf",2016) => 0.16,
         ("AQ","con",2015) => 0.17,
         ("AQ","con",2016) => 0.18,
)
[vol[r,ft,y]  = vol[r,ft,y-1] * (1-m[r,ft,y-1]) for r in regions, ft in forTypes, y in years[2:end]  ]
[println("$r \t $ft \t $y \t $(vol[r,ft,y])") for r in regions, ft in forTypes, y in years ]

LO 	 hf 	 2015 	 100.0
AQ 	 hf 	 2015 	 500.0
LO 	 con 	 2015 	 300.0
AQ 	 con 	 2015 	 700.0
LO 	 hf 	 2016 	 89.0
AQ 	 hf 	 2016 	 425.0
LO 	 con 	 2016 	 261.0
AQ 	 con 	 2016 	 581.0
8

For what it’s worth

julia> regions = ["LO","AQ"];

julia> years = collect(2015:2016);

julia> forTypes = ["hf","con"];

julia> vol = Dict(("LO","hf",2015) => 100.0,
                  ("LO","con",2015) => 300.0,
                  ("AQ","hf",2015) => 500.0,
                  ("AQ","con",2015) => 700.0,
       );

julia> m = Dict(("LO","hf",2015) => 0.11,
                ("LO","hf",2016) => 0.12,
                ("LO","con",2015) => 0.13,
                ("LO","con",2016) => 0.14,
                ("AQ","hf",2015) => 0.15,
                ("AQ","hf",2016) => 0.16,
                ("AQ","con",2015) => 0.17,
                ("AQ","con",2016) => 0.18,
       );

julia> using Base.Cartesian

julia> iter = (regions, forTypes, years[2:end]);

julia> function nloop!(vol::Dict{Tuple{I1,I2,I3},T}, m::Dict{Tuple{I1,I2,I3},T}, iter::Tuple{Vector{I1},Vector{I2},Vector{I3}}) where {I1,I2,I3,T}
         @nloops 3 i d -> begin
           iter[d]
         end begin
          vol[i_1, i_2, i_3] = vol[i_1,i_2,i_3-1] * (1-m[i_1,i_2,i_3-1])
         end
         vol
       end
nloop! (generic function with 1 method)

julia> function listcomp!(vol::Dict{Tuple{I1,I2,I3},T}, m::Dict{Tuple{I1,I2,I3},T}, regions::Vector{I1}, forTypes::Vector{I2}, years::Vector{I3}) where {I1,I2,I3,T}
         [vol[r,ft,y]  = vol[r,ft,y-1] * (1-m[r,ft,y-1]) for r in regions, ft in forTypes, y in years[2:end]  ]
         vol
       end
listcomp! (generic function with 1 method)

julia> nloop!(vol, m, iter)
Dict{Tuple{String,String,Int64},Float64} with 8 entries:
  ("AQ", "con", 2015) => 700.0
  ("LO", "hf", 2015)  => 100.0
  ("LO", "hf", 2016)  => 89.0
  ("AQ", "con", 2016) => 581.0
  ("AQ", "hf", 2016)  => 425.0
  ("LO", "con", 2015) => 300.0
  ("AQ", "hf", 2015)  => 500.0
  ("LO", "con", 2016) => 261.0

julia> #reset vol
       vol = Dict(("LO","hf",2015) => 100.0,
                  ("LO","con",2015) => 300.0,
                  ("AQ","hf",2015) => 500.0,
                  ("AQ","con",2015) => 700.0,
       );

julia> listcomp!(vol, m, regions, forTypes, years)
Dict{Tuple{String,String,Int64},Float64} with 8 entries:
  ("AQ", "con", 2015) => 700.0
  ("LO", "hf", 2015)  => 100.0
  ("LO", "hf", 2016)  => 89.0
  ("AQ", "con", 2016) => 581.0
  ("AQ", "hf", 2016)  => 425.0
  ("LO", "con", 2015) => 300.0
  ("AQ", "hf", 2015)  => 500.0
  ("LO", "con", 2016) => 261.0

julia> using BenchmarkTools

julia> @benchmark nloop!($vol, $m, $iter)
BenchmarkTools.Trial: 
  memory estimate:  1.47 KiB
  allocs estimate:  82
  --------------
  minimum time:     2.561 μs (0.00% GC)
  median time:      2.649 μs (0.00% GC)
  mean time:        2.793 μs (3.11% GC)
  maximum time:     213.791 μs (94.72% GC)
  --------------
  samples:          10000
  evals/sample:     9

julia> @benchmark listcomp!($vol, $m, $regions, $forTypes, $years)
BenchmarkTools.Trial: 
  memory estimate:  2.91 KiB
  allocs estimate:  114
  --------------
  minimum time:     2.945 μs (0.00% GC)
  median time:      3.062 μs (0.00% GC)
  mean time:        3.348 μs (6.17% GC)
  maximum time:     211.250 μs (95.48% GC)
  --------------
  samples:          10000
  evals/sample:     9

You could try benchmarking at sizes closer to what you’re using.
However, readability is definitely important too.