Create d-dimensional grid with n points in each direction (where d can be higher than 3)

New to Julia
1

I am trying to write a function that creates a d-dimensional grid (d can be higher than 3) with n points in each direction

function grid(d, n)
    return coordinates of n^d points
end

So far I have the following:

julia> using Base.Iterators

julia> n = 10
10

# 1D
julia> a = range(0, 1, length=n)
0.0:0.1111111111111111:1.0

# 2D
julia> b = collect(product(a, a))
10×10 Array{Tuple{Float64,Float64},2}:
 (0.0, 0.0)       (0.0, 0.111111)       …  (0.0, 1.0)     
 (0.111111, 0.0)  (0.111111, 0.111111)     (0.111111, 1.0)
 (0.222222, 0.0)  (0.222222, 0.111111)     (0.222222, 1.0)
 (0.333333, 0.0)  (0.333333, 0.111111)     (0.333333, 1.0)
 (0.444444, 0.0)  (0.444444, 0.111111)     (0.444444, 1.0)
 (0.555556, 0.0)  (0.555556, 0.111111)  …  (0.555556, 1.0)
 (0.666667, 0.0)  (0.666667, 0.111111)     (0.666667, 1.0)
 (0.777778, 0.0)  (0.777778, 0.111111)     (0.777778, 1.0)
 (0.888889, 0.0)  (0.888889, 0.111111)     (0.888889, 1.0)
 (1.0, 0.0)       (1.0, 0.111111)          (1.0, 1.0)     

julia> f(x) = sum(xi^2 for xi in x)
f (generic function with 1 method)

julia> f.(a)  # 1D
10-element Array{Float64,1}:
 0.0                 
 0.012345679012345678
 0.04938271604938271 
 0.1111111111111111  
 0.19753086419753085 
 0.308641975308642   
 0.4444444444444444  
 0.6049382716049383  
 0.7901234567901234  
 1.0                 

julia> f.(b)  # 2D
10×10 Array{Float64,2}:
 0.0        0.0123457  0.0493827  0.111111  …  0.604938  0.790123  1.0    
 0.0123457  0.0246914  0.0617284  0.123457     0.617284  0.802469  1.01235
 0.0493827  0.0617284  0.0987654  0.160494     0.654321  0.839506  1.04938
 0.111111   0.123457   0.160494   0.222222     0.716049  0.901235  1.11111
 0.197531   0.209877   0.246914   0.308642     0.802469  0.987654  1.19753
 0.308642   0.320988   0.358025   0.419753  …  0.91358   1.09877   1.30864
 0.444444   0.45679    0.493827   0.555556     1.04938   1.23457   1.44444
 0.604938   0.617284   0.654321   0.716049     1.20988   1.39506   1.60494
 0.790123   0.802469   0.839506   0.901235     1.39506   1.58025   1.79012
 1.0        1.01235    1.04938    1.11111      1.60494   1.79012   2.0    

julia> f.(c)  # 3D???

Is this the right approach? And if so, how can I (efficiently) extend this to higher dimensions (3 and above)?

1 Like
2

Yes, using Iterators.product is a viable approach.

That said, are you sure you want to represent grid itself directly, not eg some values at the gridpoints? The former just contains a lot of redundant information for no good reason.

3

Sometimes I will read the grid from a file but in this case I don’t necessarily need the grid itself.

Other parts of my code use Halton points:

using HaltonSequences
julia> d = 4
       p = HaltonPoint(d, length=n)
10-element HaltonPoint{Float64}:
 [0.5, 0.3333333333333333, 0.2, 0.14285714285714285]                    
 [0.25, 0.6666666666666666, 0.4, 0.2857142857142857]                    
 [0.75, 0.1111111111111111, 0.6000000000000001, 0.42857142857142855]    
 [0.125, 0.4444444444444444, 0.8, 0.5714285714285714]                   
 [0.625, 0.7777777777777777, 0.04, 0.7142857142857142]                  
 [0.375, 0.2222222222222222, 0.24000000000000002, 0.8571428571428571]   
 [0.875, 0.5555555555555556, 0.44, 0.02040816326530612]                 
 [0.0625, 0.8888888888888888, 0.6400000000000001, 0.16326530612244897]  
 [0.5625, 0.037037037037037035, 0.8400000000000001, 0.30612244897959184]
 [0.3125, 0.37037037037037035, 0.08, 0.4489795918367347]

Would it be better to store the points as a vector of vectors as with the Halton points above rather than a multidimensional array? This would be easy to read from a file but how could I generate the grid in this format?

4

However, I do need the coordinates of the points for plotting and evaluating functions at those points.

5

Would something like this help?

julia> flat_gridpoints(grids) = vec(collect(Iterators.product(grids...)))
flat_gridpoints (generic function with 1 method)

julia> a = 0:1
0:1

julia> flat_gridpoints((a, ))
2-element Array{Tuple{Int64},1}:
 (0,)
 (1,)

julia> flat_gridpoints((a, a))
4-element Array{Tuple{Int64,Int64},1}:
 (0, 0)
 (1, 0)
 (0, 1)
 (1, 1)

julia> flat_gridpoints((a, a, a))
8-element Array{Tuple{Int64,Int64,Int64},1}:
 (0, 0, 0)
 (1, 0, 0)
 (0, 1, 0)
 (1, 1, 0)
 (0, 0, 1)
 (1, 0, 1)
 (0, 1, 1)
 (1, 1, 1)

You can also collect the grid coordinates as a static vector or a vector.

3 Likes
6

I think ProductIterator this is already almost the right data type/abstraction to pass around grids in a program already, maybe together with an accessor function

factors(p::ProductIterator) = p.iterators
7

That’s perfect!

julia> creategrid(d, n) = vec(collect(Iterators.product((0:1/max(1, (n-1)):1 for _ in 1:d)...)))
creategrid (generic function with 1 method)

julia> creategrid(3,2)
8-element Array{Tuple{Float64,Float64,Float64},1}:
 (0.0, 0.0, 0.0)
 (1.0, 0.0, 0.0)
 (0.0, 1.0, 0.0)
 (1.0, 1.0, 0.0)
 (0.0, 0.0, 1.0)
 (1.0, 0.0, 1.0)
 (0.0, 1.0, 1.0)
 (1.0, 1.0, 1.0)

julia> creategrid(4,2)
16-element Array{NTuple{4,Float64},1}:
 (0.0, 0.0, 0.0, 0.0)
 (1.0, 0.0, 0.0, 0.0)
 (0.0, 1.0, 0.0, 0.0)
 (1.0, 1.0, 0.0, 0.0)
 (0.0, 0.0, 1.0, 0.0)
 (1.0, 0.0, 1.0, 0.0)
 (0.0, 1.0, 1.0, 0.0)
 (1.0, 1.0, 1.0, 0.0)
 (0.0, 0.0, 0.0, 1.0)
 (1.0, 0.0, 0.0, 1.0)
 (0.0, 1.0, 0.0, 1.0)
 (1.0, 1.0, 0.0, 1.0)
 (0.0, 0.0, 1.0, 1.0)
 (1.0, 0.0, 1.0, 1.0)
 (0.0, 1.0, 1.0, 1.0)
 (1.0, 1.0, 1.0, 1.0)
8

What is the difference between p and p.iterators and when should I use p.iterators instead of just p?

1 Like
9

p is the actual product iterator. p.iterators you would need only if you have a product p of some iterators and you want to know what these iterators are.

julia> product(1:10, 1:10).iterators
(1:10, 1:10)
2 Likes
10

Thanks for all your help. I realised I needed a vector of vectors instead of tuples so this is what I have now:

function creategrid(d::Integer, n::Integer)

    @assert d >= 1 ("d (number of dimensions) must be a positive integer")
    @assert n >= 2 ("n (number of points) must be a at least 2")

    r = range(0, 1, length = n)

    iter = Iterators.product((r for _ in 1:d)...)

    return vec([collect(i) for i in iter])
end

Any suggestions for improvement are welcome.

2 Likes