How to perform N-loops with N as an input parameter?

New to Julia
1

Hello,

Context: I would like to do a d-dimensional midlepoint integration with d as an input parameter. This leads to the problem that I need to make a d-loop

for i_1 in eachindex(x_1)
   for i_2 in eachindex(x_2)
      ...
      for i_d in eachindex(x_3)
      ...
      end
   end
...
end

which I have no idea how to perform.

The other possibility is to make an array of the d-length vectors. However, I do not know how to obtain such structure from x_1, …, x_d, either. What should be the best way to go through all n_1*…*n_d points?

Thank you in advance!

P.S. Also, I need to make it as straightforward as possible since one part of exercise is to compare the performance with MC integration for different d. Hence, no tricks to avoid loops are needed I think.

1 Like
2

This might be what you are looking for:

3

Thank you! With this, I found Iterators.jl package. There, one can find such function

product(xs...)

which returns all combinations in the cartesian product of the inputs. For instance,

x = [1,2,3]
y = [1,2,3]

for i in product(x,y)
  @show i
end

i = (1,1)
i = (2,1)
i = (3,1)
i = (1,2)
i = (2,2)
i = (3,2)
i = (1,3)
i = (2,3)
i = (3,3)

This is exactly what I need. The only problem now is that I do not know the number of input arrays in advance. Hence, I would like to add them in product one by one, and once again I am stuck. Maybe you know the solution?

Thank you in advance!

4

Isn’t CartesianRange from the linked post exactly what you need for this?

5

Not exactly I believe, or at least I do not understand how to use it the way I need. Honestly, I am a bit confused now. So, I have d-dimensional unit cube which I need to slice to cubes with edge length 1/n. Then, I need to sum over all of this cubes.

Therefore, if I want to use CartesianRange, I need to initialize d-dimensional square array A or initialize a tuple t = (n, n,..., n) of length d, but I don’t know how to do this if d is a parameter. Suppose I did this step. Now I can do something like

for i in CartesianRange(size(A))
...
end

Now, all I need is to map my indices to coordinates. Something like (k1,k2) => 1/(2*n)*(2*k1-1, 2*k2-1). If i was a tuple, it would be easy to do with collect. However, in this case i would be of CartesianIndex type, and I am once again confused.

6

To make a tuple of a given length you can use the ntuple function, e.g.

julia> ntuple(i -> [i, i+1], 10)
([1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11])

You can use this something like the following:

julia> d = 3
3

julia> myrange = CartesianRange(ntuple(x->2, d))
CartesianRange{CartesianIndex{3}}(CartesianIndex{3}((1,1,1)),CartesianIndex{3}((2,2,2)))

julia> for i in myrange
          println(i)
       end
CartesianIndex{3}((1,1,1))
CartesianIndex{3}((2,1,1))
CartesianIndex{3}((1,2,1))
CartesianIndex{3}((2,2,1))
CartesianIndex{3}((1,1,2))
CartesianIndex{3}((2,1,2))
CartesianIndex{3}((1,2,2))
CartesianIndex{3}((2,2,2))
7

Thank you very much! One last question if you don’t mind. As I mentioned, I need to make some transformation from Cartesian indices to some actual coordinates. Right now I run such routine

# d and n are input parameters

t = ntuple(i -> n, d)
myrange = CartesianRange(t)

for i in myrange
i = collect(i.I)
x = i/n - 1/(2*n)
y = f(x)
...
end

Here, f(x) is just some function of coordinates.

Is there anything one can significantly change here? I mean, it looks a bit sloppy, in my opinion.

8

You can make it a bit cleaner by not doing the collect step:

x = [(ii - 1/2) / n for ii in i.I]

Ccing @timholy in case there’s a better solution, since this would seem to be a common pattern.

9

I believe this approach is significantly slower though? I did integration routine twice: first with collect step and then without. Hence, the only difference is in the part we discussed. Here is the benchmark result for n = 6, d = 5:

BenchmarkTools.Trial: 
  memory estimate:  7.00 mb
  allocs estimate:  209955
  --------------
  minimum time:     5.754 ms (0.00% GC)
  median time:      7.358 ms (0.00% GC)
  mean time:        7.760 ms (9.68% GC)
  maximum time:     13.954 ms (22.68% GC)
  --------------
  samples:          644
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

and

BenchmarkTools.Trial: 
  memory estimate:  5.46 mb
  allocs estimate:  194403
  --------------
  minimum time:     81.565 ms (0.00% GC)
  median time:      84.981 ms (0.00% GC)
  mean time:        85.887 ms (0.69% GC)
  maximum time:     100.674 ms (0.00% GC)
  --------------
  samples:          59
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
10

If your primary concern is to perform the integration, then check out my package SparseGrids.
It uses the Iterators package for generating the D dimensional quadrature point from the 1D points.