art with code

2009-06-14

Haskell OpenGL utilities

In order to get back into the gear for writing this Haskell OpenGL application, I'll write here a small library of OpenGL utilities for roughly OpenGL 3.0 -style code. That is, no fixed-function stuff, each drawable is a struct with a shader program, streams (a.k.a. attributes), textures and uniforms.

Maybe something like the following (only works on Haskell OpenGL 2.2.3 and above.)
The code for Models and Shaders compiles but I haven't tested it. The snippets for loading textures and VBOs [see below] are working code.

First some matrix helpers (this is a snippet from a ~100-line matrix math lib):

module Matrix where
import Data.List
import Graphics.Rendering.OpenGL
import Foreign.Ptr

type Matrix4x4 = [Vec4]
type Vec4 = [GLfloat]

glMatrix :: Matrix4x4 -> IO (GLmatrix GLfloat)
glMatrix m = newMatrix ColumnMajor $ flatten m :: IO (GLmatrix GLfloat)

withMatrix4x4 :: Matrix4x4 -> (MatrixOrder -> Ptr GLfloat -> IO a) -> IO a
withMatrix4x4 matrix m = do
mat <- glMatrix matrix
withMatrix mat m

flatten :: [[a]] -> [a]
flatten = foldl1 (++)

Then the drawable definition. Drawables are structs with a program, uniforms, streams and samplers. You could think of them as curried GPU function calls, kinda like Drawable = GPU ().

module Models where
import Graphics.Rendering.OpenGL
import VBO
import Texture
import Foreign.Ptr (Ptr, castPtr)
import Matrix
import Data.Int

data Vbo a = Vbo (NumArrayIndices, BufferObject, (IntegerHandling, VertexArrayDescriptor a))
data Stream a = Stream (AttribLocation, Vbo a)
data Sampler = Sampler (UniformLocation, TextureTarget, TextureObject)

data UniformSetting = UniformSetting (UniformLocation, UniformValue)
data UniformValue =
UniformMatrix4 (Matrix4x4)
| UniformVertex4 (Vertex4 GLfloat)
| UniformVertex3 (Vertex3 GLfloat)
| UniformVertex2 (Vertex2 GLfloat)
| UniformFloat (TexCoord1 GLfloat)
| UniformInt (TexCoord1 GLint)

data Drawable = Drawable {
program :: Program,
uniforms :: [UniformSetting],
streamMode :: PrimitiveMode,
streams :: [Stream GLfloat],
indices :: Maybe (Vbo GLushort),
samplers :: [Sampler]
}

uniformSetting :: UniformSetting -> IO ()
uniformSetting (UniformSetting(location, UniformMatrix4 value)) =
withMatrix4x4 value (\order ptr -> uniformv location 16 (castPtr ptr :: Ptr (TexCoord1 GLfloat)))
uniformSetting (UniformSetting(location, UniformVertex4 value)) = uniform location $= value
uniformSetting (UniformSetting(location, UniformVertex3 value)) = uniform location $= value
uniformSetting (UniformSetting(location, UniformVertex2 value)) = uniform location $= value
uniformSetting (UniformSetting(location, UniformFloat value)) = uniform location $= value
uniformSetting (UniformSetting(location, UniformInt value)) = uniform location $= value

drawDrawable :: Drawable -> IO ()
drawDrawable d = do
currentProgram $= Just (program d)
setUniforms (uniforms d)
setSamplers (samplers d)
withStreams (streams d) (do
case indices d of
Just (Vbo (num, vbo, (_,VertexArrayDescriptor numcomp datatype stride ptr))) -> do
bindBuffer ElementArrayBuffer $= Just vbo
drawElements (streamMode d) num datatype ptr
bindBuffer ElementArrayBuffer $= Nothing
Nothing -> drawArrays (streamMode d) 0 (minNum (streams d)))
currentProgram $= Nothing
where minNum streams = minimum $ map (\(Stream (_,Vbo(n,_,_))) -> n) streams

withStreams :: [Stream a] -> IO () -> IO ()
withStreams streams m = do
setStreams streams
m
disableStreams streams

setStreams :: [Stream a] -> IO ()
setStreams streams =
mapM_ (\(Stream (location, Vbo (_, vbo, value))) -> do
bindBuffer ArrayBuffer $= Just vbo
vertexAttribArray location $= Enabled
vertexAttribPointer location $= value) streams

disableStreams :: [Stream a] -> IO ()
disableStreams streams =
mapM_ (\(Stream (location,_)) -> vertexAttribArray location $= Disabled) streams

setUniforms :: [UniformSetting] -> IO ()
setUniforms uniforms = mapM_ uniformSetting uniforms

setSamplers :: [Sampler] -> IO ()
setSamplers samplers =
mapM_ (\(i, Sampler (location, texType, tex)) -> do
activeTexture $= TextureUnit i
textureBinding texType $= Just tex
uniform location $= TexCoord1 i) $ zip [0..] samplers

Then we also need loaders for shaders, textures and buffers. Shaders are pretty easy, this loader's copied from the OpenGL binding examples:

module Shaders where
import Control.Monad
import Graphics.Rendering.OpenGL
import Graphics.UI.GLUT

loadProgram :: FilePath -> FilePath -> IO Program
loadProgram vertexShader fragmentShader =
loadProgramMulti [vertexShader] [fragmentShader]

loadProgramMulti :: [FilePath] -> [FilePath] -> IO Program
loadProgramMulti vertexShaders fragmentShaders = do
vs <- mapM readAndCompileShader vertexShaders
fs <- mapM readAndCompileShader fragmentShaders
createProgram vs fs

-- Make sure that GLSL is supported by the driver, either directly by the core
-- or via an extension.
checkGLSLSupport :: IO ()
checkGLSLSupport = do
version <- get (majorMinor glVersion)
unless (version >= (2,0)) $ do
extensions <- get glExtensions
unless ("GL_ARB_shading_language_100" `elem` extensions) $
ioError (userError "No GLSL support found.")

readAndCompileShader :: Shader s => FilePath -> IO s
readAndCompileShader filePath = do
src <- readFile filePath
[shader] <- genObjectNames 1
shaderSource shader $= [src]
compileShader shader
reportErrors
ok <- get (compileStatus shader)
infoLog <- get (shaderInfoLog shader)
mapM_ putStrLn ["Shader info log for '" ++ filePath ++ "':", infoLog, ""]
unless ok $ do
deleteObjectNames [shader]
ioError (userError "shader compilation failed")
return shader

createProgram :: [VertexShader] -> [FragmentShader] -> IO Program
createProgram vs fs = do
[program] <- genObjectNames 1
attachedShaders program $= (vs, fs)
linkProgram program
reportErrors
ok <- get (linkStatus program)
infoLog <- get (programInfoLog program)
mapM_ putStrLn ["Program info log:", infoLog, ""]
unless ok $ do
deleteObjectNames [program]
ioError (userError "linking failed")
return program

Buffers aren't much of a problem either:

module VBO where
import Foreign.Storable
import Data.Array.Storable
import Foreign.Ptr
import Graphics.Rendering.OpenGL
import Graphics.UI.GLUT

createVBO :: Storable a => [a] -> IO BufferObject
createVBO elems = do
[vbo] <- genObjectNames 1
bindBuffer ArrayBuffer $= Just vbo
arr <- newListArray (0, length elems - 1) elems -- Data.Array.MArray
withStorableArray arr (\ptr -> -- Data.Array.Storable
bufferData ArrayBuffer $= (ptrsize elems, ptr, StaticDraw))
bindBuffer ArrayBuffer $= Nothing
reportErrors
return vbo
where ptrsize [] = toEnum 0
ptrsize x:xs = toEnum $ length elems * (sizeOf x)

offset x = plusPtr nullPtr x
-- for use with e.g. VertexArrayDescriptor 3 Float 0 $ offset 0

For textures I'm using Cairo and Gdk.Pixbuf. Turning Cairo surfaces into Ptrs edible by texImage2D is a bit of a bother but eh.

module Texture where
import Data.ByteString (ByteString)
import Data.ByteString.Internal (toForeignPtr)
import Directory (doesFileExist)
import Foreign.ForeignPtr (withForeignPtr)
import Foreign.Ptr
import Graphics.Rendering.OpenGL
import Graphics.Rendering.Cairo hiding (rotate, identityMatrix)
import Graphics.UI.Gtk.Gdk.Pixbuf
import Graphics.UI.Gtk.Cairo

loadTextureFromFile :: FilePath -> IO TextureObject
loadTextureFromFile filepath = do
assertFile filepath
createTexture Texture2D Enabled
(withImageSurfaceFromPixbuf filepath $ texImage2DSurface Nothing 0)

withImageSurfaceFromPixbuf :: FilePath -> (Surface -> IO a) -> IO a
withImageSurfaceFromPixbuf filepath m = do
pixbuf <- pixbufNewFromFile filepath
w <- pixbufGetWidth pixbuf
h <- pixbufGetHeight pixbuf
withImageSurface FormatARGB32 w h (\s -> do
renderWith s (do setSourcePixbuf pixbuf 0 0
setOperator OperatorSource
paint)
m s)

assertFile :: FilePath -> IO ()
assertFile filepath = do
fex <- doesFileExist filepath
if not fex
then fail (filepath ++ " does not exist")
else return ()

createTexture :: TextureTarget -> Capability -> IO () -> IO TextureObject
createTexture target mipmap m = do
[tex] <- genObjectNames 1
textureBinding target $= Just tex
texture target $= Enabled
textureFilter target $= ((Linear', Nothing), Linear')
textureWrapMode target S $= (Repeated, Clamp)
textureWrapMode target T $= (Repeated, Clamp)
generateMipmap target $= mipmap
m
if mipmap == Enabled
then textureFilter target $= ((Linear', Just Linear'), Linear')
else return ()
textureBinding target $= Nothing
return tex

texImage2DSurface :: Maybe CubeMapTarget -> Level -> Surface -> IO ()
texImage2DSurface cubemap level imageSurface = do
pixelData <- imageSurfaceGetData imageSurface
(w,h) <- renderWith imageSurface $ do
w <- imageSurfaceGetWidth imageSurface
h <- imageSurfaceGetHeight imageSurface
return (fromIntegral w :: GLsizei, fromIntegral h :: GLsizei)
texImage2DByteString cubemap level RGBA8 w h BGRA UnsignedByte pixelData

texImage2DByteString :: Maybe CubeMapTarget
-> Level
-> PixelInternalFormat
-> GLsizei
-> GLsizei
-> PixelFormat
-> DataType
-> ByteString
-> IO ()
texImage2DByteString cubemap level iformat w h format ptype bytestring = do
let (fptr, foffset, flength) = toForeignPtr bytestring
if (fromIntegral flength) /= w * h * 4
then fail "imageSurface dimensions don't match data length"
else return ()
withForeignPtr fptr $ \ptr -> do
let optr = plusPtr ptr foffset
texImage2D cubemap NoProxy
level iformat (TextureSize2D w h) 0
(PixelData format ptype optr)

2009-06-13

Sum representation of integer powers

The power xn where x,n ∈ N can be represented as an nth-order sum where the innermost factor is n!.

First, consider x1: the difference between x1 and (x+1)1 is 1. So we have a first-order sum where the innermost factor is 1! (= 1.)

x^1 = sum(i = 0..x, 1)

Now let's look at x2. Let's write the sequence of x2 from -3 to +3 and mark under each two elements the difference between them:

9 4 1 0 1 4 9
-5 -3 -1 1 3 5
2 2 2 2 2

From this we see that we have a second-order sum where the innermost factor is 2! (=2) and the outer factor is -1 (going diagonally left from zero.)

x^2 = -1*sum(i = 0..x, 2*i*sum(j = 0..i, 1))

We can now give a recursive function for the sums that takes a list of factors as its argument (this one only works for positive integers though):

sumMap lst f = sum (map f lst)

sums [] i = 0
sums [x] i = x * i
sums (x:xs) i = x*i + sumMap [0..i] (sums xs)

Then, let's look at the left side of x6:

46656 15625 4096 729 64 1 0 1
-31031 -11529 -3367 -665 -63 -1 1
19502 8162 2702 602 62 2
-11340 -5460 -2100 -540 -60
5880 3360 1560 480
-2520 -1800 -1080
720 720

From which we get the factors [-1, 62, -540, 1560, -1800, 720] and

x^6 = -1*sum(a=0..x, 62*sum(b=0..a, -540*sum(c=0..b, 1560*sum(d=0..c, -1800*sum(e=0..d, 720*sum(f=0..e, 1))))))
or
pow6 = sums [-1, 62, -540, 1560, -1800, 720]

In the general case, the first factor is -1(n-1), the second factor is -2n + 2 * -1n, the second to last factor is -(n-1)n! * 2-1 and the last factor is n!. I don't know about the rest of the factors.

Fun with cubes


Here's the difference grid for the third power:

-27 -8 -1 0 1 8 27
19 7 1 1 7 19
-12 -6 0 6 12
6 6 6 6

Which gives the factors [1, -6, 6].

One result of the sum representation is that each x cubed minus x is divisible by 6: 6 | x3 - x. Or put another way, x3 = x + 6k, k ∈ Z. And the sum between two cubes minus the sum of the bases is also divisible by six: 6 | x3 + y3 - (x + y).

We also see that the difference between any two integer cubes grows as a second order function: n3 - (n-1)3 = 6(n/2)(n+1)+1 = 3n2+3n+1.

cubeDiff n = 3*(n^2 + n) + 1
pow3 n = sumMap [0..n-1] cubeDiff

-- pow3 5
-- 125

There's another amusing difference grid in (x3 - x) / 6:

0 1 4 10 20 35 56 84
1 3 6 10 15 21 28
2 3 4 5 6 7
1 1 1 1 1

2009-06-10

A wee bit of algebra

Am currently in the process of writing a game of algebra in Haskell. The gameplay consists of coaxing group axioms to equality or, in the case of the neutral element and the inverse function, a binding. The base act of equation munging is accomplished through one's big bag of previously proven lemmas and given axioms, which one then applies to matching parts of the equation in the feeble hope of achieving enlightenment.

For example, let us step through the process of proving associativity for a o b = (a x b) x 1, given the abelian group x. The | x in the following denotes applying the rule x at the emboldened position in the equation (why yes, it is an AST internally):

a o (b o c) = (a o b) o c | a o b = (a x b) x 1
(a x (b o c)) x 1 = (a o b) o c | a o b = (a x b) x 1
(a x ((b x c) x 1)) x 1 = (a o b) o c | (a x b) x c = a x (b x c)
(a x (b x (c x 1))) x 1 = (a o b) o c | (a x b) x c = a x (b x c)
((a x b) x (c x 1)) x 1 = (a o b) o c | a x b = b x a
((a x b) x (1 x c)) x 1 = (a o b) o c | (a x b) x c = a x (b x c)
(((a x b) x 1) x c) x 1 = (a o b) o c | a o b = (a x b) x 1
((a o b) x c) x 1 = (a o b) o c | a o b = (a x b) x 1
(a o b) o c = (a o b) o c

Oh, have the neutral element as well:

a o e(o) = a | a o b = a x b x 1
(a x e(o)) x 1 = a | (a = b) = ((a x inv(x,1)) = (b x inv(x,1)))
((a x e(o)) x 1) x inv(x,1) = a x inv(x,1) | (a x b) x c = a x (b x c)
(a x e(o)) x (1 x inv(x,1)) = a x inv(x,1) | a x inv(x,a) = e(x)
(a x e(o)) x e(x) = a x inv(x,1) | a x e(x) = a
a x e(o) = a x inv(x,1) | (a = b) = ((inv(x,c) x a) = (inv(x,c) x a))
inv(x,a) x (a x e(o)) = inv(x,a) x (a x inv(x,1)) | (a x b) x c = a x (b x c)
(inv(x,a) x a) x e(o) = inv(x,a) x (a x inv(x,1)) | inv(x,a) x a = e(x)
e(x) x e(o) = inv(x,a) x (a x inv(x,1)) | (a x b) x c = a x (b x c)
e(x) x e(o) = (inv(x,a) x a) x inv(x,1) | inv(x,a) x a = e(x)
e(x) x e(o) = e(x) x inv(x,1) | e(x) x a = a
e(x) x e(o) = inv(x,1) | e(x) x a = a
e(o) = inv(x,1)

As you might imagine, it is not a very exciting game in its current form. One might imagine replacing the typographical symbols with mushrooms and small furry animals, and the addition of fireworks and showers of colourful jewels on the successful completion of one's proof obligations. Maybe even a fountain of money.

Roam the countryside and befriend local wildlife with your awesome powers of axiomatic rigor! Bring down the castle drawbridge with a well-placed induction! Banish ghosts with the radiant aura of reductio ad absurdum!

Reducing the rest of mathematics to gameplay mechanics remains an open question.

Codewise, what I'm doing looks like this (leaving out all the hairy bits):

type Op = String
data Expr = Expr (Op, Expr, Expr)
| A | B | C
| Inv (Op, Expr)
| Neutral Op
| Literal Int

type Pattern = Expr
type Substitution = Expr

data Rule = Rule (Pattern, Substitution)

type CheckableRule = ((Rule -> Bool), Rule)


isTrue :: Rule -> Bool
isTrue (Rule (a,b)) = a == b

isBinding :: Expr -> Rule -> Bool
isBinding e (Rule (a,b)) = e == a || e == b

opf :: Op -> Expr -> Expr -> Expr
opf op a b = Expr (op, a, b)


{-
Group axioms main stage:

Stability: a o b exists in G for all a, b in G
Magma!

Associativity: a o (b o c) = (a o b) o c
Semigroup!

Neutral element: a o 0 = 0 o a = a
Monoid!

Inverse element: a o a_inv = a_inv o a = 0
Group! Enter Abelian bonus stage!
-}
type CheckableRule = ((Rule -> Bool), Rule)

checkCheckableRule :: CheckableRule -> Bool
checkCheckableRule (p, rule) = p rule

associativity :: Op -> CheckableRule
associativity op = (isTrue, (A `o` (B `o` C)) `eq` ((A `o` B) `o` C))
where o = opf op

rightNeutral :: Op -> CheckableRule
rightNeutral op = (isBinding (Neutral op), (A `o` Neutral op) `eq` A)
where o = opf op

leftNeutral :: Op -> CheckableRule
leftNeutral op = (isBinding (Neutral op), (Neutral op `o` A) `eq` A)
where o = opf op

rightInverse :: Op -> CheckableRule
rightInverse op = (isBinding (Inv (op, A)), (A `o` Inv (op, A)) `eq` Neutral op)
where o = opf op

leftInverse :: Op -> CheckableRule
leftInverse op = (isBinding (Inv (op, A)), (Inv (op, A) `o` A) `eq` Neutral op)
where o = opf op

{-
Abelian bonus stage:

Commutativity: a o b = b o a
Abelian group! Enter ring bonus stage!
-}

commutativity :: Op -> CheckableRule
commutativity op = (isTrue, (A `o` B) `eq` (B `o` A))
where o = opf op

magma op = [] -- FIXME?
semiGroup op = magma op ++ [associativity op]
monoid op = semiGroup op ++ [rightNeutral op, leftNeutral op]
group op = monoid op ++ [rightInverse op, leftInverse op]
abelianGroup op = group op ++ [commutativity op]

{-
Ring bonus stage:

Show bonus function to be a semigroup!

Distributivity of x over o:
a x (b o c) = (a x b) o (a x c)
(a o b) x c = (a x c) o (b x c)
Pseudo-ring!
-}

leftDistributivity :: Op -> Op -> CheckableRule
leftDistributivity opO opX = (isTrue, (A `x` (B `o` C)) `eq` ((A `x` B) `o` (B `x` C)))
where o = opf opO
x = opf opX

rightDistributivity :: Op -> Op -> CheckableRule
rightDistributivity opO opX = (isTrue, (A `x` (B `o` C)) `eq` ((A `x` B) `o` (B `x` C)))
where o = opf opO
x = opf opX
{-
Neutral element for x: a x 1 = 1 x a = a
Ring!

Commutativity for x: a x b = b x a
Commutative ring! Enter field bonus stage!

Field bonus stage:

Inverse element for x in G: a x a_inv = a_inv x a = 1
Field! Superior! Shower of jewels!
-}

pseudoRing o x = abelianGroup o ++ semiGroup x ++
[leftDistributivity o x, rightDistributivity o x]
ring o x = pseudoRing o x ++ [rightNeutral x, leftNeutral x]
commutativeRing o x = ring o x ++ [commutativity x]
field o x = commutativeRing o x ++ [rightInverse x, leftInverse x]

2009-06-05

A moment of extrapolation

Hitachi exec tells that by 2014, spinning rust will have 15TB of capacity (and I guess around 140MB/s bandwidth and 7ms seek times.)

Extrapolating the flash SSD curve with a yearly doubling in density, you'd get around 1TB of flash for 100e in 2014. And not only has flash density doubled yearly, the bandwidth has done the same.

Current flash SSDs boast around 200MB/s of bandwidth. Assuming it continues to double every year, the 1TB flash drive of 2014 would have 6.4GB/s bandwidth.

I don't know if latency is keeping pace, but assuming the 75us latencies of today fall 20% annually, they'd be at 25us. But who knows.

Anyhow, 6.4GB/s IO bandwidth in a cheapo computer would be awesome.

Blog Archive