Struggling with a Point type for several dimensions

New to Julia
21

You can use Setfield.jl and do p2 = @set p1.x += 1 which is equivalent to (i.e., as efficient as) p2 = Point3D(p1.x + 1, p1.y, p1.z).

22

You can use StructArrays.jl for SOA vs AOS:

array = StructArray{Point{2, Float32}}(rand(Float32, 100), rand(Float32, 100))
# and then just use `array` wherever you'd use `Vector{Point{2, Float32}}` and it will magically all work with SOA :)
23

Thanks again

24

Just wanted to point out that immutable 3d/2d points being the most performant isnā€™t unique to Julia. Iā€™ve written plenty of C++ physics/geometry processing code and Iā€™ve always treated my point structs as immutable. Even for SoA you often want to write to a new buffer and then simply swap pointers after ( for random access pattern, serial update can be ok for mutation ). Cache canā€™t perform as well if it knows itā€™s data might be invalidated by a write.

25

@Tamas_Papp kind of mentioned it briefly, but you may also want to check out in-place operators, usually denoted with a !

It might look something like

function updatePoint!(p1::Array{Point3D, 1}, i, nx, ny, nz)
  p1[i] = Point3d(p1[i].x + nx, p1[i].y + ny, p1[i].z + nz)
  end;

for updating a single point i. Then to move to broadcasting

function updatePoints!(p1::Point3D, nx, ny, nz)
  p1 = Point3d(p1.x + nx, p1.y + ny, p1.z + nz)
  end;

and call as updatePoints!.( vector/matrix of point3d, same sized nx, ny, nz ) to update many at once

Add: Ah, points .= updatePoints!.(points, ā€¦ in this case?

26

How do the other approaches perform?

struct Point3D ##<: FieldVector{3, Float64}
    x::Float64
    y::Float64
    z::Float64
end

pv = Point3D.(collect(1.0:0.1:100000),1.0,1.0)
n = length(pv)
xp,yp,zp = rand(n), rand(n), rand(n)
@time for iter=1:74 
  pv .= updatePoints!.(pv, xp,yp,zp);
  end;

Gives 0.99sec for near 74,000,000 simple point updates

27

Btw, you can also write this roughly as:

pv = rand(SVector{3, Float64}, n)
pv .= pv .+ SVector.(xp, yp, zp)
1 Like
28

Just a smidge slower, ~0.993 to ~1.003 sec with the FieldVector part uncommented and using SVector

using StaticArrays

sv = SVector.(xp, yp, zp)
@time for iter=1:74 pv .= pv .+ sv; end;

And even closer with the pv = line, ~0.997

I think @time was pretty accurate, i did use a few runs

BenchmarkTools.Parameters(5.0, 10000, 1, 0.0, true, false, 0.05, 0.01)
times: [1.00113e9, 1.00184e9, 1.0025e9, 1.00361e9, 1.00452e9]
mem: 3552, alloc: 148
BenchmarkTools.Parameters(5.0, 10000, 1, 0.0, true, false, 0.05, 0.01)
times: [9.89299e8, 9.8971e8, 9.90509e8, 9.91875e8, 9.91978e8, 9.92714e8]
mem: 4736, alloc: 148
BenchmarkTools.Parameters(5.0, 10000, 1, 0.0, true, false, 0.05, 0.01)
times: [1.00491e9, 1.00499e9, 1.00532e9, 1.00674e9, 1.00691e9]
mem: 3552, alloc: 148

with loopsize 1
x = @benchmark for iter=1:1 pv .= pv .+ sv; end;
times: ā€œmin: 1.3489198e7, mean: 1.3560357623306233e7ā€ (13.x millis)
mem: 48, alloc: 2, samples: 369

x2 = @benchmark for iter=1:1 pv .= px.updatePoints!.(pv, xp,yp,zp); end;
times: ā€œmin: 1.3358128e7, mean: 1.3417005728494624e7ā€
mem: 64, alloc: 2, samples: 372

x3 = @benchmark for iter=1:1 pv .= pv .+ sv; end; -Fieldvec
times: ā€œmin: 1.3453972e7, mean: 1.3529067089189189e7ā€
mem: 48, alloc: 2, samples 370

29

You should probably use @btime from Benchmarktools to time that, since itā€™s not in a function and you are using a global variable.

1 Like
30

Putting it in one line is pretty important here, since that will then fuse the calculation together, without allocating a temporary Vector{Point} :wink:
Also, thatā€™s how I would benchmark things;


using StaticArrays
pv = SVector.(collect(1.0:0.1:100000),1.0,1.0)
n = length(pv)
xp,yp,zp = rand(n), rand(n), rand(n)
function update_point(p1::SVector, nx, ny, nz)
    SVector(p1[1] + nx, p1[2] + ny, p1[3] + nz)
end

# move the expression you want to benchmark into function scope 
# (optional with benchmark tools and $ infront of input params)
test1(pv, xp, yp, zp) = pv .= pv .+ SVector.(xp, yp, zp)
test2(pv, xp, yp, zp) = pv .= update_point.(pv, xp, yp, zp)

using BenchmarkTools
# If you just want the time
@btime test1($pv, $xp, $yp, $zp)
@btime test2($pv, $xp, $yp, $zp)

# if you want to compare results + and have it statistically sound!
b1 = @benchmark test1($pv, $xp, $yp, $zp)
b2 = @benchmark test2($pv, $xp, $yp, $zp)

judge(minimum(b1), minimum(b2)).time == :invariant

Btw, I think StaticArrays should deprecate fieldvectorsā€¦ Now, that one can overload getfield (getproperty), one can also get svector.x without a struct.
If anyone would like to have that, Iā€™d be open to accept a PR to GeometryTypes adding getproperty overloads for the Point type (which is just a StaticVector inheriting form StaticArrays) :wink:

2 Likes
31

include("points2.jl")

13.316 ms (0 alloc, 0 byte) 13.332 ms minimum? in b1
13.321 ms (0 alloc, 0 byte) 13.322 ms minimum for b2
true