# How to efficiently skip calculating some vector elements?

There are significant efficiency gains by implementing Projective Geometric Algebra (PGA) in OpenGL shader hardware as described in Steven De Keninck’s recent paper, “Look, Ma, No Matrices!

Less significant (because it applies to only a few applications), I’d like to implement a more efficient Julia port of his C++ Projective Geometric Algebra reference implementation by efficiently skipping the calculation of some vector elements that are not needed in some applications. For example, the following code detects the intersection of two line segments, using only the first element of the calculated vector. In this application, calculating the other elements (i.e., 2:8) is a waste.

``````# mwe.jl - Minimal Working Example
# detect intersecting line segments

# define multivector basis names
# 0 denotes projective dimension (e.g., e0 * e0 = 0)
basis = [	# iField
"1"		#  1	scalar
"e0"	#  2	grade 1 vectors
"e1"	#  3
"e2"	#  4
"e01"	#  5	grade 2 vectors (bivectors)
"e20"	#  6
"e12"	#  7
"e012"]	#  8	pseudoscalar

# define the basis elements
nField = 2^3 # 3 = 2D + 1 dimensions
e0 =   zeros(Float32, nField); e0[2] = 1
e1 =   zeros(Float32, nField); e1[3] = 1
e2 =   zeros(Float32, nField); e2[4] = 1
e01 =  zeros(Float32, nField); e01[5] = 1
e20 =  zeros(Float32, nField); e20[6] = 1
e12 =  zeros(Float32, nField); e12[7] = 1
e012 = zeros(Float32, nField); e012[8] = 1

# regressive product: vee operator (&, \vee)
function Base.:&(a::Vector{Float32},b::Vector{Float32})::Vector{Float32}
res = similar(a)
res[1]=a[1]*b[8]+a[2]*b[7]+a[3]*b[6]+a[4]*b[5]+a[5]*b[4]+a[6]*b[3]+a[7]*b[2]+a[8]*b[1]
res[2]=a[2]*b[8]+a[8]*b[2]+a[5]*b[6]-a[6]*b[5]
res[3]=a[3]*b[8]+a[8]*b[3]-a[5]*b[7]+a[7]*b[5]
res[4]=a[4]*b[8]+a[8]*b[4]+a[6]*b[7]-a[7]*b[6]
res[5]=a[5]*b[8]+a[8]*b[5]
res[6]=a[6]*b[8]+a[8]*b[6]
res[7]=a[7]*b[8]+a[8]*b[7]
res[8]=a[8]*b[8]
return res
end

function point(x::Number, y::Number)::Vector{Float32}
return x*e20 + y*e01 + e12
end

function tst()
for iTest = 1:2
P0A = point(0,0)
P0B = point(2,0)
P1A = point(1,1)
P1B = iTest == 1 ?
point(0,-1) :
point(-2,-1)

L0 = P0A & P0B # Line 0
A0A = L0 & P1A # oriented area about Line 0
A0B = L0 & P1B # oriented area about Line 0
L1 = P1A & P1B # Line 1
A1A = L1 & P0A # oriented area about Line 1
A1B = L1 & P0B # oriented area about Line 1

println("iTest: \$iTest")
println("L0:  \$L0")
println("A0A: \$A0A")
println("A0B: \$A0B")
println("L1:  \$L1")
println("A1A: \$A1A")
println("A1B: \$A1B")

isIntersection =
(A0A[1] * A0B[1] < 0) &&
(A1A[1] * A1B[1] < 0) ? true : false
println("isIntersection: \$isIntersection\n")
end # iTest
end # tst()
``````
``````julia> include("mwe.jl")
tst (generic function with 1 method)

julia> tst()
iTest: 1
L0:  Float32[0.0, 0.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0]
A0A: Float32[-2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A0B: Float32[2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
L1:  Float32[0.0, 1.0, -2.0, 1.0, 0.0, 0.0, 0.0, 0.0]
A1A: Float32[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A1B: Float32[-3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
isIntersection: true

iTest: 2
L0:  Float32[0.0, 0.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0]
A0A: Float32[-2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A0B: Float32[2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
L1:  Float32[0.0, -1.0, -2.0, 3.0, 0.0, 0.0, 0.0, 0.0]
A1A: Float32[-1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A1B: Float32[-5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
isIntersection: false
``````

In C or C++, I would reverse the order of vector element calculations and use a switch statement (without break lines) so that a single test would jump to the code to calculate the first n (e.g., n=1 in the line segment intersection code above) elements. However, Julia doesn’t have switch statements and it doesn’t seem that Julia’s macros were intended to “ungenerate” code. Is there an efficient way to skip the calculation of some vector elements?

If you use StaticArrays.jl, which is probably a good idea anyway since the dimensions of your vectors are probably fixed for a given set of calculations, then there’s a good chance that dead code elimination by the compiler will strip out the computations of elements you don’t subsequently use.

That being said, there are various packages that give you a switch-like statement in Julia, though it will generally compile down to a sequence of `if` branches at the end.

1 Like

Thanks for the suggestions. Currently the vector dimensions are fixed, but that is because I haven’t yet been able to port ganja.js’s ability to dynamically change vector dimensions as described in this demonstration of the FABRIK inverse kinematics algorithm at SIGGRAPH 2019. However, my port is close enough to working that I suspect Julia can also dynamically change PGA vector dimensions.

Maybe this package is what you’re looking for: GitHub - JuliaAPlavin/Skipper.jl

They don’t have to be fixed for all time to use StaticArrays effectively. They just need to be determined early in the computation, not changing in your inner loops.