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

GDP projections

I took some GDP and population figures from Wikipedia and NationMaster and tried to predict the future with them.

First off, the current situation:

Region2008 nominal GDP/capitaPopulation in millions
United States47,103310
European Union38,390500
Japan38,055127
Russia12,487142
China3,1741,345
India1,0781,173


Grabbing last decade's average growth rates for each, we can do a linear extrapolation to 2020:

Region2020 GDP/capita2020 population
United States72,900341
European Union73,200661
Japan46,000127
Russia91,000137
China13,3811,431
India2,9901,368


These projections do have some issues though:

Russia posted 6x growth over the last decade, continuing that for a second decade sounds unrealistic, because 1998 Russia was suffering from a big financial crisis. A more realistic 3-4x decade growth would put them at $37,500-$50,000 GDP per capita.

160 million person more people in the EU. The EU population growth is mostly driven by enlargement. To grow by 160 million, the EU would need the Balkans (25 million), Turkey (75 million) and Ukraine (46 million.) And the EUR/USD exchange rate is favoring EUR at the moment, while in 2000 the USD was higher. Using an estimate based on 1999 figures, the GDP/capita projection would be around $67000. Though if the EU population growth misses the spot, the GDP/capita will likely be higher as a result, the possible member countries having a GDP/capita a third or fourth of the EU average.

Japan's low growth is from the bubble bursting in 1996, and the following decade of flat growth. Japan's only now reaching mid-90s levels again; the 1994 and 2008 Japan GDP/capita are roughly the same. Assuming that Japan starts posting decade growth rates in the 1.5-1.75x range, they'd be somewhere between the US and EU in 2020.

2040 projections


The linear projections fail even more for 2040. To highlight, consider Russia's GDP/capita. The 2040 projection is 3.7 million dollars. The whole Russian GDP would be 478 trillion dollars, constituting half the global economy. Which is unlikely, to say the least.

My guess for Russia's growth over the next three decades would be something like 4x in 2010-2020, 2x in 2020-2030, 1.5x in 2030-2040. Which'd put Russia at $156,000 GDP/capita. It might seem high today, but consider that the US GDP/capita thirty years ago was $10,000.

At the rate EU and Indian populations were growing in the last decade, EU's population would overtake India in 2080, and they'd have 3.5 and 3.4 billion people, respectively. The 2040 projection for EU would be 1.15 billion citizens, and that'd require all of North Africa and Middle East, plus perhaps parts of Sub-Saharan Africa, and Russia or Pakistan. Which sounds a wee bit far-fetched. Perhaps there will be an union of unions, or some other export mechanism for the EU's rule-of-law-investment-free-trade-combo.

The EU GDP/capita in 2040 would be around $153,000. The US GDP/capita projection for 2040 is $165,000, and Japan might be somewhere around the two. The linear projection for Japan would have them at $78,000, but that doesn't sound too likely.

The Chinese GDP/capita would be at $174,262 at today's high growth rates, which I wouldn't trust to hold for three decades. Today China's GDP/capita compared to the US is at around 1970s South Korea, or late-50s Japan and Spain. If China sticks to Japan's path (4x, 4.5x, 2.7x), they'll be at around US levels in 2040, perhaps followed by a bubble crash. With the South Korean way (4x, 2.3x, 1.6x), the projection is 35% of the US GDP/capita, or $58,000. The Spanish path (2x, 5x, 2x) would take them to half the US GDP/capita.

India would be the most populous country in the world with 1.86 billion, propelling it past China's projected 1.6 billion. India's 3x per decade growth would bring their GDP/capita to around $18,800. The Indian growth might pick up though, three decades of Japan-style growth would bring them to $58,000. Perhaps a good estimate would be somewhere in between: $32,000?

And then the computational power disclaimer: by 2040, you can get a human brain worth of processing power for something like $100 / year, which nicely throws a spanner in these GDP projections. If each person has a dozen virtual humans running on the side, the amount of stuff that one person can get done might well be multiplied by that.

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]

Blog Archive