Archive for category Snippets

Snippet: Class that never raises AttributeError

Sometimes we’d like to have a class that is tolerant of attribute accesses that weren’t anticipated at class design time, returning a default ‘None’ value for an uninitialized attribute. Using defaultdict to override the `__dict__` attribute makes this very easy:

from collections import defaultdict

class X(object):
    def __init__(self):
        self.__dict__ = defaultdict(lambda : None)
    def __getattr__(self, attr):
        return self.__dict__[attr]

Now we can initialize an object with some attributes, but access others and get back the default None:

x = X()
x.color = "RED"
x.size = 6

print x.color
print x.size
print x.material

prints:

RED
6
None

Of course, this defeats a key validation feature in Python, so you’ll need to take extra care that you specify attributes correctly.

Leave a comment

Snippet: Uniquify a sequence, preserving order

I came across some code today that used a set to keep track of previously seen values while iterating over a sequence, keeping just the not-seen-before items. A brute force kind of thing:

seen = set()
unique = []
for item in sequence:
    if not item in seen:
        unique.append(item)
        seen.add(item)

I remembered that I had come up with a simple form of this a while ago, using a list comprehension to do this in a single expression. I dug up the code, and wrapped it up in a nice little method. This version accepts any sequence or generator, and makes an effort to return a value of the same type as the input sequence:

def unique(seq):
    """Function to keep only the unique values supplied in a given 
       sequence, preserving original order."""
       
    # determine what type of return sequence to construct
    if isinstance(seq, (list,tuple)):
        returnType = type(seq)
    elif isinstance(seq, basestring):
        returnType = type(seq)('').join
    else:
        # - generators and their ilk should just return a list
        returnType = list

    try:
        seen = set()
        return returnType(item for item in seq if not (item in seen or seen.add(item)))
    except TypeError:
        # sequence items are not of a hashable type, can't use a set for uniqueness
        seen = []
        return returnType(item for item in seq if not (item in seen or seen.append(item)))

My first pass at this tried to compare the benefit of using a list vs. a set for seen – it turns out, both versions are useful, in case the items in the incoming sequence aren’t hashable, in which case the only option for seen is a list.

Leave a comment

RK integrator

While moving files from my old laptop drive to my new one, I found a nice Runge-Kutta integrator class that I had written ages ago. So long ago, in fact, that I was a little embarrassed at the newbiness of some of my code. So I decided to update my code to get a nice RK class out of it, using list comprehensions instead of “for i in range” loops, and including an integrate method that acts as a generator so that the calling code can cycle through each integration step. As is typical in R-K, the system state is maintained in a vector X, and the calling method must provide a callback function that will return dX/dt.

Here is the class:

class RKIntegrator :
    "Class used to perform Runge-Kutta integration of set of ODE's"

    def __init__( self, dt, derivFunc, degree=0, initConds=None ):
        self.dt = float(dt)
        self.dt_2 = dt / 2.0
        self.t = float(0)
        if not (degree or initConds):
            raise ValueError("must specify degree or initial conditions")
        if initConds is not None:
            self.x = initConds[:]
        else:
            self.x = [0.0 for i in range(degree)]
        self.derivFunc = derivFunc


    def doIntegrationStep( self ):
        dt = self.dt
        dxFunc = self.derivFunc
        t2 = self.t + self.dt_2

        dx = dxFunc( self.t, self.x )
        delx0 = [ dx_i*dt for dx_i in dx ]
        xv = [x_i + delx0_i/2.0 for x_i, delx0_i in zip(self.x, delx0)]
        dx = dxFunc( t2, xv )
        delx1 = [ dx_i*dt for dx_i in dx ]
        xv = [x_i + delx1_i/2.0 for x_i,delx1_i in zip(self.x,delx1)]
        dx = dxFunc( t2, xv )
        delx2 = [ dx_i*dt for dx_i in dx ]
        xv = [x_i + delx1_2 for x_i,delx1_2 in zip(self.x,delx2)]
          
        self.t += dt
        dx = dxFunc(self.t, xv)
        
        self.x = [
            x_i + ( delx0_i + dx_i*dt + 2.0*(delx1_i + delx2_i) ) / 6.0
                for x_i, dx_i, delx0_i, delx1_i, delx2_i in 
                    zip(self.x, dx, delx0, delx1, delx2)
            ]
    
    def integrate(self):
        while True:
            self.doIntegrationStep()
            yield self.t, self.x

Here is an example of finding X with constant acceleration of 4:

def getDX( t, x ):
  return [ 
    x[1], 
    4.0 
    ]

isWhole = lambda x : abs(x-round(x)) < 1e6 

rk = RKIntegrator( dt=0.1, derivFunc=getDX, initConds = [0.0, 0.0] )
for t,x in rk.integrate():
    if t > 10: 
        break
    if isWhole(t):
        print t,', '.join('%.2f' % x_i for x_i in x)

Googling for ‘Python runge kutta’, I came across this blog posting:
http://doswa.com/blog/2009/01/02/fourth-order-runge-kutta-numerical-integration/
This does a good job, but hardcodes the vector size to just x, velocity, and acceleration. Here is how my R-K integrator would implement Doswa’s code:

def accel(t,x):
    stiffness = 1
    damping = -0.005
    x,v = x
    return -stiffness*x - damping*v

def getDX(t,x):
    return [
        x[1],
        accel(t,x)
        ]

rk = RKIntegrator( dt=1.0/40.0, derivFunc=getDX, initConds = [50.0, 5.0] )
for t,x in rk.integrate():
    if t > 100.1: 
        break
    if isWhole(t):
        print t,', '.join('%.2f' % x_i for x_i in x)

My results match the posted results to 2 places.

1 Comment

Code Golf – Converting seconds to hh:mm:ss.mmm, or reduce() ad absurdam

I’ve had to do some timing using time.clock(), but I like seeing output as hh:mm:ss.000.  divmod is the perfect tool for this, as it does the division and remainder in a single call (rekindling my love for Python’s ease in returning multiple values from a function):

def to_hhmmss(sec):
    min,sec = divmod(sec,60)
    hr,min = divmod(min,60)
    return "%d:%02d:%06.3f" % (hr,min,sec)

Right about the time I wrote that, I happened to see an example using reduce, and thought about using successive calls to divmod to build up the tuple to be passed in to this function’s string interpolation.  That is, I needed to call reduce in such a way to convert 3601.001 to the tuple (1, 0, 1, 1) (1 hour, 0 minutes, 1 second, and 1 millisecond).

reduce() applies a binary operation to a sequence by taking the first two items of the sequence, performing the operation on them and saving the result to a temporary accumulator, then applies the binary operation to the accumulated value and the 3rd item in the sequence, the to the 4th item, etc.  For instance to use reduce to sum up a list of integers [1,5,10,50], we would call reduce with the binary operation:

fn = lambda a,b: a+b

and reduce would work through the list as if running this explicit code:

acc = 1
acc = fn(acc, 5)
acc = fn(acc, 10)
acc = fn(acc, 50)
return acc

I need reduce to build up a tuple, which is a perfectly good thing to do with tuple addition.
Also, I’m going to convert the milliseconds field into its own integer field also, so I’ll need to convert our initial value containing seconds to a value containing milliseconds (also allowing me to round off any decimal portion smaller than 1 msec). And I’ll use divmod with a succession of divisors to get the correct time value for each field.

So the sequence of values that I will pass to reduce will be the succession of divisors for each field in the tuple, working right to left.  To convert to hh, mm, ss, and msec values, these divisors are (1000, 60, 60), and these will be the 2nd thru nth values of the input sequence – the 1st value will be the value in milliseconds to be converted.

The last thing to do is to define our binary function, that will perform the successive divmods, will accumulating our growing tuple of time fields.  Maybe it would be easiest to just map out how our sample conversion of 3601.001 seconds (or 3601001 msec) would work, with some yet-to-be-defined function F.  The last step in the reduce would give us:

(1,0,1,1) = F((0,1,1), 60)
  (0,1,1) = F((1,1), 60)
    (1,1) = F(msec, 1000)

Well, that last line (representing the first binary reduction) looks bad, since all the others are taking a tuple as the first value.  It seems that each successive reduction grows that value by one more term, so the first term should be a 1-value tuple, containing our value in milliseconds:

    (1,1) = F((msec,), 1000)

Also, it’s not really true that this step would have to return (1,1).  Just so long as the final tuple ends with (…,1,1).  So I’ll redo the sequence of steps showing X’s for the unknown intermediate values:

(1,0,1,1) = F((X,1,1), 60)
  (X,1,1) = F((X,1), 60)
    (X,1) = F((msec,), 1000)

The heart of our mystery function F is essentially a call to divmod, using the 0’th element in the current accumulator:

F = lambda a,b : divmod(a[0],b)

But F has to also carry forward (accumulate) the previous divmod remainders, found in the 1-thru-n values in the tuple in a.  So we modify F to include these:

F = lambda a,b : divmod(a[0],b) + a[1:]

Now we have just about all the pieces we need. We have our the binary function F.  The sequence of values to pass to F is composed of the initial accumulator (msec,), followed by the divisors (1000, 60, 60).  We just need to build our call to reduce, and use it to interpolate into a suitable hh:mm:ss-style string:

def to_hhmmss(s):
    return "%d:%02d:%02d.%03d" % \
        reduce(lambda a,b : divmod(a[0],b) + a[1:], [(s*1000,),1000,60,60])

This may not be the most easy-to-read version of to_hhmmss, but it is good to stretch our thinking into some alternative problem-solving approaches once in a while.

, ,

Leave a comment