Interactive graphics - Python, SymPy and VisPy

Submitted by rybu on Wed, 02/10/2016 - 00:11

One of my childhood jobs was for an engineering company.  This was back in the days before Matlab or Mathematica.  I would write software for this company (interestingly the company failed and was reincarnated as FlexPipe in 2001, the owners being the previous engineers), to compute numerical integrals to predict the pressure it should take to make their pipe explode, in various configurations.   Eventually I used the Persistence of Vision raytracer to create slick graphics depicting the internal structure of their pipe, to help sell the product to potential investors. 

Raytracing was extremely time consuming in the 1990's so I developed a fairly efficient pipeline for generating the graphics.  One that I've been using ever since.   Below is an image created by this pipeline, about 5 years ago. My pipeline primarily involved writing a Pascal (eventually C++) script to generate (part of) the script for Persistence of Vision.  I would then run Persistence of Vision, sometimes for a few days, to generate the final image.

 

The main advantage of using Pascal and C++ to generate the data is it allowed me to be highly prescriptive on the geometric objects used -- allowing for the computation of the precise location of the circles in the above picture, as these are circles that the intersect the knot in precisely five points.  The primary disadvantage is writing the Pascal or C++ program is the relatively slow compilation process and transition to a final rendered object.  I've been sitting on the fence for years wondering if I should use an interactive modelling environment like AutoCad or Rhino. 

I believe I've stumbled across a potential shortcut that will help me avoid that.

Python is a popular high-level programming language that has a large variety of increasingly sophisticated libraries available.  VisPy is a library that utilizes OpenGL to interface with a computers GPU. In short, this allows for the rapid production of software to create rapidly-rendered graphics.  So rapid, in many cases, that interactive graphics are possible.   

It would be nice if it was that straightforward.  I've been coding in Python to see how possible all this is.   Using a high-level language like Python comes with a price.  So I've set myself a set of exercises to see if I can code "the way I'd like to" in a high-level way, and get Python to do what I like.  This is the first in a series of posts where I hope to convince myself that I can. 

The first obstacle I've encountered is the sympy library. This is Python's primary library for manipulating formal mathematical expressions, such as functions, roots, special numbers such as π, etc. It is ultimately an object that stores mathematical constructions in a tree.  So if you have a complicated mathematical expression and wish to evaluate it, sympy will parse through that tree (a structure based on pointers) and slowly assemble the required binary operations, perhaps using automatic casts and such to decide what you really want the library to do in the process.  For a decent-sized formula (like say the boundary of a tubular neighbourhood of a torus knot in R3) this can easily be tens of thousands of times slower than coding it up in C++.  If this was the end of the story, I would discard Python right here. Luckily sympy has libraries that allow one to speed this process up. 

Here is my first coding example, using vispy to generate an interactive visualization of a (p,q)-torus knot.  In future posts I hope to make the graphics increasingly interactive, using more features of the vispy library. ​

The code is broken into four blocks.

  1. Is the high-level block where everything is written as abstract sympy algebraic expressions.
  2. Is where we experiment with the various methods for converting sympy expressions into callable functions.  In this case, ufunctify appears to be the most efficient.
  3. We generate the mesh data.
  4. We pass the data to VisPy for rendering.

One other thing to note, I am running this code in an iPython / Jupyter notebook.  This provides a web-browser based environment for editing and running Python code.   I've been teaching a course recently using this environment and have found it quite pleasant.   I will put the code in a comment below. 

And here is a screenshot of the output.

In my next post I will explore the VisPy canvas a little further, and hopefully bring further layers of interactivity into this program.

 

Comments

rybu

Thu, 02/11/2016 - 13:15

A somewhat enticing option is the Theano library.  It allows users to ask the GPU to do (mostly) linear-algebra tasks, with floating points and integers.  

I tried calling Theano to generate the mesh for my (p,q) torus knots but the CPU using ufuncify was faster. 

On top of that, I found setting up Theano to work properly was a bit of a hassle.  It does not install cleanly via apt-get or pip, as one has to set up a .theanorc file.  My main problem was not the lack of documentation (I could not find any with the install) but the sheer number of competing voices when I went on-line to look it up. 

rybu

Mon, 02/29/2016 - 17:51 

## BLOCK 1
import sympy as sp

## symbols we need to describe p,q-torus knots
## t time parameter. p,q indexes the torus knot
## r minor radius, R major radius
spt, spp, spq, spr, spR = sp.symbols("t p q r R", real=True)

c = sp.Matrix([(spR+spr*sp.cos(2*sp.pi*spq*spt))*sp.cos(2*sp.pi*spp*spt), 
     (spR+spr*sp.cos(2*sp.pi*spq*spt))*sp.sin(2*sp.pi*spp*spt), 
      spr*sp.sin(2*sp.pi*spq*spt)])
dc = sp.Matrix([sp.diff(x,spt) for x in c]) # derivative
ldc = sp.sqrt(sum( [ x**2 for x in dc ] )).simplify() # speed
udc = dc/ldc

## 2nd order
kc = sp.Matrix([sp.diff(x,spt) for x in udc]) # curvature vector
ks = sp.sqrt(sum( [ x**2 for x in kc])) # curvature scalar
ukc = kc/ks # unit curvature vector

## bi-normal
bnc = udc.cross(ukc) # cross of unit tangent and unit curvature.

## the parametrization of the boundary of the width w tubular neighbourhood

spw, spu = sp.symbols("w, u", real=True) ## width of torus knot, and meridional parameter
tSurf = c + spw*sp.cos(2*sp.pi*spu)*ukc + spw*sp.sin(2*sp.pi*spu)*bnc

## BLOCK 2

## Let's have a visualization routine that takes as input a curve and a framing of the curve.  
##  We will then plot things like a tubular neighbourhood of the curve, together with some
##  decoration on the boundary. 

import numpy as np
import itertools as it

## (a) lambdify with numpy.  This returns a 3-element list.
knotSnp = sp.lambdify((spt, spp, spq, spr, spR, spw, spu), tSurf, "numpy" )

## (b) ufuncify
from sympy.utilities.autowrap import ufuncify
knotSuf = [ufuncify([spt, spp, spq, spr, spR, spw, spu], tSurf[i]) for i in range(3)]

## (c) theano
from sympy.printing.theanocode import theano_function
knotSth = theano_function([spt,spp,spq,spr,spR,spw,spu], [tSurf],
                          dims={spt:0, spp:0, spq:0, spr:0, spR:0, spw:0, spu:0})

kp = 5 ## these are the "p" and
kq = 3 ## "q" of our (p,q) torus knot.
tR = 1.6 # major torus radius
tr = 0.6 # minor torus radius. 
kt = (np.pi*tr) / (4*kp) # knot radial thickness 2*pi*tr is circumf, and kp strands pass through so this
## should be around 2*pi*tr  would be 2*kp*kt for the knot to fill the surface, i.e kt = pi*tr / 4*kp
## make bigger or smaller depending on how much empty space one wants to see.

seg = kp*300 ## segments along length of pq torus knot. kp*120 gives a fairly smooth image.
segm = 40 ## meridional segmentation of pq torus knot. 60 is fairly smooth. 

def surf1(i,j): ## sympy raw
    return np.array(tSurf.evalf(subs={spt:float(i)/seg, spu:float(j)/segm, 
                                       spp:kp, spq:kq, spr:tr, spR:tR, spw:kt}) )
def surf2(i,j): ## lambdify
    return np.array(knotSnp(float(i)/seg, kp, kq, tr, tR, kt, float(j)/segm)).ravel()
def surf3(i,j): ## ufuncify
    return np.array([knotSuf[k](float(i)/seg, kp, kq, tr, tR, kt, float(j)/segm) for k in range(3)])
def surf4(i,j): ## theano
    return knotSth(float(i)/seg, kp, kq, tr, tR, kt, float(j)/segm).ravel()

Surf = [surf1, surf2, surf3, surf4]
SurfLabel = ["sympy.evalf", "sympy.lambdify", "ufuncify", "theano"]

## BLOCK 3

k = 2 # determines which method we use to cast sympy expressions to a callable function.

import time as ti
start=ti.time()
surf = Surf[k]
pqMesh = np.array(list( it.chain.from_iterable(([surf(i,j), surf(i+1,j), surf(i,j+1)], 
                                                [surf(i+1,j), surf(i+1,j+1), surf(i,j+1)]) for i,j in it.product(range(seg), range(segm)))))
end=ti.time()
print(SurfLabel[k]+" mesh generation: "+str(end-start)+" seconds.", flush=True)

colA = np.array([0.0, 0.0, 1.0, 1.0])
colB = np.array([1.0, 1.0, 0.0, 1.0])

#print(SurfLabel[useFunc] + str(" mesh ")+str(pqMesh))

def iCol(i,j): ## interpolates between colA and colB
    alph = np.pi*float(j)/segm
    bet = np.pi*float(i)/seg
    gam = 1.0*alph - kp*kq*bet ## coprime combination of alpha and beta
    return (np.sin(gam)**2)*colA + (np.cos(gam)**2)*colB
    
## coloration of mesh
pqColors = np.array(list(it.chain.from_iterable(([iCol(i,j),iCol(i+1,j),iCol(i,j+1)], 
                                              [iCol(i+1,j),iCol(i+1,j+1),iCol(i,j+1)]) for i,j in it.product(range(seg), range(segm)) )))

## BLOCK 4

## Let's now generate the vispy data
from vispy import geometry, app, scene
import vispy, sys
from vispy.scene import cameras

canvas = scene.SceneCanvas(keys='interactive', size=(800, 600), show=True)
view = canvas.central_widget.add_view()
view.parent = canvas.scene
view.camera = scene.TurntableCamera(parent=view.scene, fov=45)  
view.bgcolor = '#606060'
## shading options are empty (pass nothing), shading='smooth' and shading='flat'
pqtorus = scene.visuals.Mesh(pqMesh,face_colors=pqColors, shading='smooth')

axis = scene.visuals.XYZAxis(parent=view.scene)

## and showtime! 
view.add(pqtorus)
canvas.app.run()