Archive for March, 2008

A defection

So my employer decided to buy me a Macbook Pro. Stay tuned for a lengthy post in which I summarise the reasonably smooth but lengthy process of setting up gfortran, scipy, matplotlib and all my other tools.

My home computers still run Ubuntu, but I think it’s safe to say there will be fewer posts here about how to do this or that in Ubuntu.


Leave a Comment

The return of Pyticles

It’s been a while, and I’ve written some code and made some decisions. The first is to make use of python’s object oriented features , and the second is not to let efficiency concerns get in the way of clean design.

In the structure I’ve chosen, a particle object contains the positions, momenta and other properties of n particles (there is no individual particle object). Forces are also objects, and are added to particle systems – this is to allow dynamic adding of multiple forces. Each force has it’s own neighbour list – this might seem inefficient until you notice that long range forces have a much larger list (that needs updating much less often) than short range forces.

Leave a Comment

Generic nbody ODE solver

So, I wrote a generic nbody ODE solver, in pure python. It’s likely quite slow, but it works – seemed like too much work at the time to interface to scipy’s ODE solvers, which seemed designed for smaller, less coupled systems.

I may have misjudged, though which is why I’m noting this link ,which has an example of what seems to be a large dimensional ODE system. Thing is, we still need to do things like update adjacency lists and distances, and recompute summation densities between steps. Worth a look, though.

Leave a Comment

Calling IDL from python

I’m reading at the moment. Yes, my goal in life is to do everything via python scripts.

Leave a Comment

random numpy wisdom – record arrays

> Is it possible to create a numpy array from an array of a C structure like
> this?
> struct RateInfo
> {
> unsigned int ctm;
> double open;
> double low;
> double high;
> double close;
> double vol;
> };

Sure. On the numpy side, you would make an record array with the
appropriate dtype and size.

In [1]: from numpy import *

In [2]: dt = dtype([(‘ctm’, uint), (‘open’, double), (‘low’, double),
(‘high’, double), (‘close’, double), (‘vol’, double)])

In [3]: a = empty(10, dtype=dt)

On the C side, you would iterate through your C array and your numpy
array and just assign elements from the one to the other. If you have
a contiguous C array, you could also just use memcpy().

This is probably reliable because all of your struct members take up
multiples of 4 bytes and most C compilers will pack those without any
space between them. If you were mixing, say, chars and doubles, the C
compiler may try to align the doubles on a 4-byte boundary (or
possibly another boundary, but 4-bytes is common). In that case, you
will have to figure out how your C compiler is packing the member and
emulate that in your dtype. Each of the tuples in the constructor can
have a third element which represents the byte offset of that member
from the beginning of the struct.

In [4]: dt2 = dtype([(‘ctm’, uint, 0), (‘open’, double, 4), (‘low’,
double, 12), (‘high’, double, 20), (‘close’, double, 28), (‘vol’,
double, 36)])

Robert Kern

Leave a Comment

IDL’s coolest feature (debuggin python)

My favourite thing about IDL is that if it encounters an error, it drops you in the middle of the program in the interpreter, so you can examine variables, try commands, and generally work your way out of the error. I’ve seen the following tip elsewhere before, but I’m posting it here mostly for my own reference. Thanks to infinite8s on the pyglet mailing list for this one:

in the main key_press function:

def on_key_press(self, symbol, modifiers):
if symbol == key.Q:
elif symbol == key.D:
import pdb

now if you run your game and hit d, you’ll be in dropped into the
debugger right after the call to set_trace() (so basically in the
middle of your event loop). From here you can usually walk the call
graph, and investigate things if need be, or set breakpoints, etc. If
you hit ‘c’ at the debugger prompt your program will continue running.
Of course, you probably want to remove this code when you finally
release your product :)

Comments (1)

Python object method versus external function speed

On a plane from Qld I came up with all sorts of clever OO ways to refactor my particle system. Unfortunately it looks as though the cost is too high. Not by a great deal, but for a realtime app that I want to squeeze as much out as possible it’s (just barealy) not acceptable.

In short, calling traverse(neighbour_list) is around 8% faster than calling neighbour_list.traverse(). It appears that it’s better (faster, anyway) to write mutator methods external to the class they need to work on.

Still thinking about this one – it’s quite possible that once all the maths involved with each interaction or force contribution is included, the fractional cost will get lower are lower (I’m prepared to accept ~%2, I think. Yeah, 2 percent is a nice, arbitrary small number).

Starting to come to my senses. 10% performance hit is actually quite acceptable, especially given this program is all about playing and prototyping. Being able to dynamically add and remove forces, and for each force to have its own independent (or shared) neighbour list is important.

Better to have a convenient design, and use techniques such as pointed out here for speeding up. The most relevent seem to be: using local variables (e.g. assign force.apply() to apply(), before the loop body) as python accesses locals much faster than globals, using map and list comprehensions where possible, and use xrange for ranges.

""" The purpose of this script is to explore the cost of expressing
things in an object oriented way.
from time import time
import Numeric
import math
n = 20
class Nlist:
def __init__(self):
self.nip = 0
self.iap = Numeric.zeros((n*n,2))
def build(self):
for i in range(n):
for j in range(i+1,n):
self.iap[self.nip,0] = i
self.iap[self.nip,1] = j
self.nip += 1
def traverse(self):
for k in range(nl.nip):
i = nl.iap[k,0]
j = nl.iap[k,1]
class Force:
def __init__(self,nl): = nl
def pairwise(self,i,j):
g = math.sqrt(i*j)
def apply(self):
for k in range(nl.nip):
i = nl.iap[k,0]
j = nl.iap[k,1]
def traverse(nl):
for k in range(nl.nip):
i = nl.iap[k,0]
j = nl.iap[k,1]
def apply_force(f,nl):
for k in range(nl.nip):
i = nl.iap[k,0]
j = nl.iap[k,1]
nl = Nlist()
print "Timing externel function for nl traversal versus object method"
start_time = time()
end_time = time()
t1 = end_time - start_time
print "object method: %f elapsed" %(t1)
start_time = time()
end_time = time()
t2 = end_time - start_time
print "external function: %f elapsed" %(t2)
diff = (t1-t2)/t1
print "Diff as fraction of t1 %f" %(diff)
print "Timing force object with (same) attached neighbour lists versus seperate neighbour list"
bigforce = Force(nl)
smallforce = Force(nl)
start_time = time()
end_time = time()
t1 = end_time - start_time
print "internal force loop: %f elapsed" %(t1)
start_time = time()
end_time = time()
t2 = end_time - start_time
print "external force loop: %f elapsed" %(t2)
diff = (t1-t2)/t1
print "Diff as fraction of t1 %f" %(diff)

Leave a Comment