The best way to define two or more arrays, similar to a given one

General Usage
1

Is there a method to avoid repeating similar(u) in the code below?
To be more precise, I define a parameterized surface (u,v)->(f(u,v), g(u,v), h(u,v)) with parameters in a rectangle [a, b] x [c,d]:

ru = a:0.1:b
rv = c:0.1:d
u, v = [ui for ui in ru, vi in rv], [vi for ui in ru, vi in rv]

x, y, z = similar(u), similar(u), similar(u)
@. x = f(u,v)
@. y = g(u,v)
@. z =h(u,v)
2

You can define something, if this happens enough:

julia> similars_v2(x...) = (similar(x...) for _ in 1:999);

julia> a,b,c = similars_v2([1,2], Float64);

julia> a[1] = 99;

julia> b
2-element Vector{Float64}:
 2.1863129326e-314
 2.2054957863e-314

Edit – initially I suggested this, which is a bad idea. Since Iterators.repeated is an ordinary function, its argument is evaluated once, and thus you get the same array repeatedly, not different arrays:

julia> similars(x...) = Iterators.repeated(similar(x...));

julia> a,b,c = similars([1,2], Float64);

julia> a[1] = 99;

julia> b
2-element Vector{Float64}:
 99.0
  2.175639858e-314
1 Like
3

I’m coming back because I’ve noticed some strange behavior in IJulia.

If I paste these lines of code

u = [1, 2, 3.0, -4, 1, 0]
v = [ 5, 1, 0, 3, 7, 2.0]

x, y, z = similars(u)
@. x = u*u - v 
@. y = 2*u+u*v
@. z = v+u-u^2

in the same cell, and inspect x, y, z, they contain the same elements.

Initially (when I accepted this solution) I wrote in a cell:

@. x=...

and inspected its elements, then in another cell:

@. y=

etc and they had distinct elements. Why when the lines are written in the same cell x, y, z seem to have the same reference and as a consequence the same values?

4

You need no meshgrids or no similars. You just need to do the following:

u = range(a, b; length=m+1)
v = range(c, d; length=n+1)

x = @. f(u, v')
y = @. g(u, v')
z = @. h(u, v')

Working example:

using Plots
pyplot()

function halftorus(n=25, a=5, b=10; lims=(-11, 11), size=(500, 400))
    u = range(-π/2, π/2, length=n+1)
    v = range(0, 2Ď€, length=2n+1)
    x = @. (b + a * cos(u')) * cos(v)
    y = @. (b + a * cos(u')) * sin(v)
    z = @.      a * sin(u')  * one(v)
    surface(x, y, z; lims, colorbar=false, size)
end

halftorus()

image

1 Like
5

Thank you @genkuroki,
I know this definition, but I need the meshgrid, because I define a triangulation on the surface. First I construct the Delaunay triangulation on the points of the meshgrid, and then I lift it to the surface through the parameterization function.

6

Oh you are right, my answer is too quick. It calls similar once, not N times. Sorry! This is exactly like the fill(ones(2), 3) mistake.

7

If u and v are meshgrids, you can make x, y, and z with no similars, more easily than above, by

x = @. f(u, v)
y = @. g(u, v)
z = @. h(u, v)

Note the positions of the @. macros.

1 Like
8

No, it doesn’t work. Commenting out:

x, y, z =similars(u)

in the enclosed lines of code or in my original code with u, v as meshgrids, generates the error:

UndefVarError: x not defined

Well, I’ll keep the less inspired setting:

x, y, z = similar(u), similar(u), similar(u)

because it works. I asked this question just for finding out a more specialized method to allocate memory simultaneously, for x, y, z.

9

Working example with no similars:

using Plots
pyplot()

ru, rv = -1.5:0.1:1.5, -1:0.1:1
u, v = [ui for ui in ru, vi in rv], [vi for ui in ru, vi in rv]

x = @. u*u - v 
y = @. 2*u+u*v
z = @. v+u-u^2

surface(x, y, z; colorbar=false, size=(500, 400))

image

Please pay attention to the positions of @.s.

Postscript: (for further clarity)

For an AbstractVector u, the code @. x = f(u) is essentially equivalent to

for i in eachindex(u)
    x[i] = f(u[i])
end

If x is not defined in this code, then we get “UndefVarError: x not defined”.

But the code x = @. f(u) is essentially equivalent to

x = [f(u[i]) for i in eachindex(u)]

We obtain the new vector x.

If we just notice this difference, we can remove the redundant similars.

Working example:

u = range(0, 1; length=7)
@. x = sinpi(u)

UndefVarError: x not defined
...

x = @. sinpi(u)

7-element Vector{Float64}:
 0.0
 0.5
 0.8660254037844386
 1.0
 0.8660254037844387
 0.4999999999999999
 0.0
2 Likes
11

Thank you so much for the solution and especially for explaning each definition!

12

Not that this is much more inspiring, but this is a possibility if that appears in another context:

x, y, z = ntuple(i->similar(u),3)

and of course this could be made into a simple function like:

similars(u,N) = ntuple(i->similar(u),N)

x, y, z = similars(u,3)
2 Likes