Symbolic Optimisation and JIT Compilation of Geometric Algebra in Python with GAALOP and Numba
Hugo Hadfield, Dietmar Hildenbrand, Alex Arsenovic
2019
GARAMON
Versor
Gaigen
GAALOP
Clifford (Python)
Matlab GA Toolbox
Ganja.js
And many others!
Low level implementations are fast to run but slow to write
High level implementations are fast to write, fun to work with, but typically slow to run
Lots of pieces in place for easy collaboration and prototyping
# Import Conformal GA (4,1)
from clifford.g3c import *
# Import prebuilt tools for conformal GA
from clifford.tools.g3c import *
# Import pyganja for visualisation
from pyganja import *
object_array = [
random_conformal_point(),
random_point_pair(),
random_circle(),
random_line(),
random_sphere(),
random_plane()
]
# The simplest way to use pyganja
draw(object_array,
static=False ,
scale=0.1)from clifford.tools.g3c import *
from clifford.tools.g3c.rotor_parameterisation import *
from pyganja import *
# Create two objects
X1 = random_circle()
X2 = random_circle()
# Find the rotor between them
R = rotor_between_objects(X1, X2)
# Take the log of the rotor
logR = general_logarithm(R)
# Blend from 0 to 1
rotor_list = [general_exp(0.1*alpha*logR)
for alpha in range(10)]
# Apply the interpolated rotors
object_list = [R_int*X1*~R_int
for R_int in rotor_list]
# The simplest way to use pyganja
draw(object_list ,
static=False ,
scale=0.1)
def example_one(X, R):
Y = 1 + (R * X * ~R)
Ystar = Y * (e1 ^ e2 ^ e3 ^ einf ^ e0)
F = Ystar | (e1 ^ e2)
return FWe would like to accelerate this, compile it perhaps, to avoid run-time type checking overhead. Ideally with no code changes:
@make_this_faster
def example_one(X, R):
Y = 1 + (R * X * ~R)
Ystar = Y * (e1 ^ e2 ^ e3 ^ einf ^ e0)
F = Ystar | (e1 ^ e2)
return F
import numba
import random
@numba.njit
def monte_carlo_pi(nsamples):
acc = 0
for i in range(nsamples):
x = random.random()
y = random.random()
if (x ** 2 + y ** 2) < 1.0:
acc += 1
return 4.0 * acc / nsamplesTime for 100,000,000
Pure Python: 76.67s
Numba JITed: 1.87s
41x faster
@numba.njit # THIS DOES NOT WORK!
def example_one(X, R):
Y = 1 + (R * X * ~R)
Ystar = Y * (e1 ^ e2 ^ e3 ^ einf ^ e0)
F = Ystar | (e1 ^ e2)
return F
Have to hand optimise to get fast code
Class support is not yet good enough for rich multivector classes
@numba.njit # Efficient but ugly
def example_one_faster_jit(X, R):
Y = gp(gp(R, X), rev(R))
Y[0] += 1
Ystar = gp(Y, e12345_value)
F = ip(Ystar, e12_value)
return F# MAYBE WE CAN DO SOMETHING LIKE THIS ?
@GAALOP
def example_one(X, R):
Y = 1 + (R * X * ~R);
Ystar = Y * (e1 ^ e2 ^ e3 ^ einf ^ e0);
F = Ystar.(e1 ^ e2);
?F;
GAALOPScript has a pretty similar syntax to Clifford Python, lets just use that!
Would also let us use existing GAALOPScript code in Python without having to rewrite the code :)
THE TRANSLITERATED CODE WE GET BACK FROM GAALOP IS JITTABLE WITH NUMBA!
from gajit import *
from clifford.tools.g3c import random_conformal_point, random_sphere
# The GAALOPScript function goes here
body = """
R = *(A ^ C ^ einf);
PP = R ^ (*S);
hasIntersection = PP.PP;
"""
# Create a python function from that
# GAALOPScript function
raytrace = Algorithm(
inputs = ['S', 'A', 'C'],
outputs = ['hasIntersection'],
intermediates = ['R','PP'],
body=body,
function_name="raytracer",
blade_mask_list = [layout.grade_mask(4),
layout.grade_mask(1),
layout.grade_mask(1)]
)
# We can now call the algorithm
# from python, passing in python
# multivectors and getting a python
# multivector back
S = random_sphere()
A = random_conformal_point()
C = random_conformal_point()
result = raytrace(S,A,C)@gajit([layout.grade_mask(4),layout.grade_mask(1),layout.grade_mask(1)])
def clifford_raytrace_gajit(S,A,C):
PP = ((e12345 * (A ^ C ^ einf)) ^ (e12345 * S))
hasIntersection = (PP | PP)
return hasIntersectionGAALOPScript has a pretty similar syntax to Clifford Python, maybe we can just translate it automatically?
We can pass grade masks on the inputs to ensure only the minimum number of symbols are generated
from clifford.g3c import *
from clifford.tools.g3c import *
from gajit import *
@gajit([layout.grade_mask(4),layout.grade_mask(1),layout.grade_mask(1)])
def clifford_raytrace_gajit(S,A,C):
PP = ((e12345 * (A ^ C ^ einf)) ^ (e12345 * S))
hasIntersection = (PP | PP)
return hasIntersection
S = random_sphere()
A = random_conformal_point()
C = random_conformal_point()
hit = clifford_raytrace_gajit(S,A,C)Not only does this allow us to use GAALOP as a precompiler for Python, it also makes Python a valid front end for GAALOP
from clifford.g3c import *
from clifford.tools.g3c import *
from gajit import *
@gajit([layout.grade_mask(4),layout.grade_mask(1),layout.grade_mask(1)],
galgebra=True)
def clifford_raytrace_gajit(S,A,C):
PP = ((e12345 * (A ^ C ^ einf)) ^ (e12345 * S))
hasIntersection = (PP | PP)
return hasIntersection
S = random_sphere()
A = random_conformal_point()
C = random_conformal_point()
hit = clifford_raytrace_gajit(S,A,C)@gajit([layout.grade_mask(4),layout.grade_mask(1),layout.grade_mask(1)])
def clifford_raytrace_gajit(S,A,C):
PP = ((e12345 * (A ^ C ^ einf)) ^ (e12345 * S))
hasIntersection = (PP | PP)
return hasIntersectionRequires no changes to source code. Generates Highly efficient implementations.
Can target multiple output languages/frameworks.