# Which approach for multidimensional computation?

**URL:** <https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182>\
**Category:** Data\
**Tags:** question, data\
**Created:** [August 2, 2017, 3:23pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182 "2017-08-02T15:23:19Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [August 2, 2017, 3:23pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/1 "2017-08-02T15:23:19Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 2, 2017, 3:36pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/2 "2017-08-02T15:36:43Z")

</div>

> [@sylvaticus](#):
>
> 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.

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

> **[GitHub - ahwillia/Einsum.jl: Einstein summation notation in Julia](https://github.com/ahwillia/Einsum.jl)**
>
> Einstein summation notation in Julia. Contribute to ahwillia/Einsum.jl development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [August 2, 2017, 4:35pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/3 "2017-08-02T16:35:30Z")

</div>

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

```julia
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)

```

```julia
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 

```

---

<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:** [August 2, 2017, 6:40pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/4 "2017-08-02T18:40:14Z")

</div>

Cartesian.jl is flexible: [https://github.com/timholy/Cartesian.jl](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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 2, 2017, 6:47pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/5 "2017-08-02T18:47:42Z")

</div>

> [@Elrod](#):
>
> Cartesian.jl is flexible: [GitHub - timholy/Cartesian.jl: Fast multidimensional algorithms](https://github.com/timholy/Cartesian.jl)

The modern version of that is Base.Cartesian

[https://docs.julialang.org/en/latest/devdocs/cartesian/](https://docs.julialang.org/en/latest/devdocs/cartesian/)

---

<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:** [August 2, 2017, 6:52pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/6 "2017-08-02T18:52:23Z")

</div>

Ah, oops. Using Base.Cartesian:

```julia
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
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.

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [August 3, 2017, 9:40am UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/7 "2017-08-03T09:40:31Z")

</div>

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

```julia
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

```

---

<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:** [August 3, 2017, 2:47pm UTC](https://discourse.julialang.org/t/which-approach-for-multidimensional-computation/5182/8 "2017-08-03T14:47:21Z")

</div>

For what it’s worth

```julia
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.