Debugging C/API python extensions

I feel really fracking dumb right about now. But at least I now know how to use GDB (the Gnu DeBugger) on this kind of project. So I've got that going for me. Which is nice.

In which we discover that gcc is not Latin

Had a very odd bug today while installing Prover9. Figured I'd record it for my (and maybe others') future reference. (The current version as of this writing is 2011-11A; at least one older version had the same behavior on my machine.)

Background: I'm messing around currently on a lightweight Linux called Puppy. Like most (all?) flavors of Linux, one is encouraged, when adding software, to compile it from source when possible. I've used gcc, the out-of-the-box GNU compiler, to compile a couple of little C++ exercises; but not until today had I decided to compile anything nontrivial.

Prover9 is designed to be downloaded as a tarball, unpacked, and then compiled without requiring any smarts on the part of the user. The unpacking was easy.

The process of compilation was supposed to be accomplished by simply giving the command

make all

to the shell. Now, keep in mind that at the start of this process I have only the slightest idea what is supposed to happen when I issue this command.

Much to my frustration, after a big bunch of output lines had scrolled across my terminal, make halted and caught fire. The topmost error in the stack, as much headscratching finally elucidated, was

undefined reference to round

while gcc was trying to compile one of the C sources.

Google sent me to this Stackoverflow question, which was thankfully not voted down (though a few people apparently tried). I still don't quite understand the details: basically, older versions of C core libraries don't come with certain obvious functions (in this case, round and ceil), and getting the compiler to pull in appropriate definitions for these functions requires issuing some extra instructions.

What's interesting is that the author of the compile (make) script issued (one version of) these extra instructions. If you unpack the archive and open up /path/to/archive/LADR-2009-11A/provers.src/Makefile, you will see at line 66

prover9: prover9.o $(OBJECTS)
    $(CC) $(CFLAGS) -lm -o prover9 prover9.o $(OBJECTS) ../ladr/libladr.a

The part that I'm interested in is the -lm instruction flag. It basically is needed in the course of linking libraries (like the math library containing round). However, my version of gcc wants this flag to be at the end of the line -- it complains (and more importantly, compilation halts and catches fire) with the flag in its current location.

The user can of course go into a text editor and manually edit Makefile so that -lm follows $(OBJECTS) (separated by a space on both sides). This worked for me: compilation completed successfully and all three included tests passed.

Triangulations...

... or how I learned that eventually I'll have to stop worrying and love a compiler.

On objects

I haven't had cause to blog about it much at all, but the only course I'm taking this semester is a for-fun course on approximation theory, or more specifically splines.

The guy teaching it, Larry, is... shall we say, old-school. I don't mean that in any way negatively; just that he remembers the days of punch-cards first-hand. However, I do mean to say that his way of computer programming is very different from my own -- at least, as I am now. The computational part of the course has given me an opportunity to reflect on my own changing thoughts, intuitions and preferences when it comes to programming.

Subclassing immutable types in Python 3

I've been hacking around in Python 3 for a while now, writing (as I mentioned a while back) a package for implementing arbitrary finite first-order structures.

The natural way to build such a thing is to subclass sets in some way: the traditional way of describing a first-order structure is as a set with additional information attached. But Python has two native set types: set and frozenset, the former mutable and the later im-.

Ideally, we'd like to enforce that you can't add or subtract elements from a model -- mean, if you plan to iterate through a model M, it would be a poor choice to pop elements off one at a time and throw them away, so a good programmer should make that an action that a user can't do by accident. Hence, I'm going to subclass from frozenset instead of the mutable set class.

There's only one problem: initializing an immutable object, or a mutable object which inherits from an immutable class, requires doing things a little differently.

I won't go through the whole setup process for my Model class, since that would involve a lot of explaining, but I'll do a simplified example: a derived class from frozenset with an extra data attribute foo.

But first, where does the problem crop up anyway? Let's say that we didn't know about this whole business: how would we usually program a class inheritance?

class SetWithFoo(frozenset):
    def __init__(self,X,foo_in):
        frozenset.__init__(self,X)
        self.foo = foo_in

 But when you compile this code we might get something like

>>> S = SetWithFoo({1,2},"bar")
>>> S.foo
    'bar'
>>>1 in S
    False

What's gone wrong? Well, remember how frozensets are immutable? And remember how __init__(self,...) isn't a constructor, because the object self already exists? What that means here is that self gets summoned into existence from the void with certain elements -- and those are the only elements which will ever belong to it. By the time __init__ sees the object, it can't change its members.

The method which does the summoning is where we need to work. That method is called __new__(...), and it's the only method in your class which doesn't take self as an argument -- because self doesn't exist yet! Instead, the first argument to __new__ is the class of the object it's creating. You the programmer don't have to worry about this at all -- just put cls as the first argname, and Python will take care of the rest:

def __new__(cls,X,foo_in):

Now, at this point there are two schools of thought on what to do next. One of these schools says that if you're going to bother overriding __new__, you should code the whole initialization in there and just leave __init__ alone (don't even explicitly override it). I'm more in the other side, which thinks that only that which has to be done in __new__ (that is, what has to be done before freezing the basic data of your object, in this case the immutable members of the set) should be done there; everything else can be handled profitably in __init__. The one caveat is that the arglists (including default arguments) of the two methods must be the same (except for cls and self, of course), or else Python will throw a fit when you call SetWithFoo(args)
Long story short, either of the two following code blocks will do just fine.

def __new__(cls,X,foo_in):
    s = frozenset.__new__(cls,X)
    s.foo = foo_in

or

def __new__(cls,X,foo_in):
    return frozenset.__new__(cls,X)

def __init__(self,X,foo_in)
    self.foo = foo_in

but obviously not both ;)

Blogging learning Python: Equivalence Relations II

Part I

Last time we looked at the basic class methods of the eqrel class. Python knows these are class methods because their definition blocks are indented under the main class header. What this means in practice is that a method like union is called on an instance EQ of the equivalence relation class, modifies that instance, and doesn't return anything.

The methods we'll see this time are different: they don't belong to the class but to the overall module, and they return a new instance of the class.

As mentioned, the two main functions we want to be able to compute are the equivalence relation join and meet. Join doesn't require any new technology:

def eqjoin(EQa,EQb):
    """Returns an instance of eqrel coding the least equivalence relation
    containing both EQa and EQb.

    Both arguments must be of the same length."""
    if len(EQa) != len(EQb): raise IndexMismatchError()
    
    N = len(EQa)
    EQr = eqrel(N)

    for iter1 in range(N):
        for iter2 in [y for y in [EQa[iter1],EQb[iter1]] if y >= 0]:
            EQr.union(iter1,iter2)

    return EQr

I've included a custom error message in case __name__ passes stupid data to the method; I'm sure that there are other possible errors I could anticipate, but this is the only one I can remember accidentally tripping over in practice. Here's the custom exception definition:

class IndexMismatchError(Exception):
    def __init__(self):
        self.msg = "Index Mismatch!\n\nAll equivalence relations must have the same length!"

    def __str__(self):
        return repr(self.msg)

Notice that a Python exception is a class (which for obvious reasons inherits from Exception). Basically all I've written this one to do is print a message detailing why it got called; nothing fancy.

Now let's program equivalence relation meet. But to do this, we need a low-wattage implementation of breadth-first search. (Recall that we already have a class method EQ.alists(), shown in the last post, for creating an adjacency-list representation of the tree out of the original parent data.)

def bfs_component(alists,origin):
    """bfs_component(alists) returns a list, the connected component
    of the element origin in the graph coded by alists."""
    N = len(alists)
    is_searched = [False] * N
    queue = [origin]
    ret_component = []

    while len(queue) > 0:
        active_element = queue.pop()
        if is_searched[active_element]:
            pass
        else:
            is_searched[active_element] = True
            ret_component.append(active_element)

            for iter1 in [x for x in alists[active_element] if not is_searched[x]]:
                queue.insert(0,iter1)

    return ret_component

All this implementation does is compute the connected component of whatever vertex number is passed as the origin argument. This will allow us to search downwards, where the array-of-parents data representation allows easy upward search but no downward search. (The tradeoff is that it is more difficult to tend a tree represented as adjacency lists, cut off and reattach elements, etc. Not an insurmountable problem, but annoying.)

Finally we can meet the meet:

def eqmeet(EQa,EQb):
    """eqmeet(EQa,EQb) returns an instance of eqrel coding the greatest
    equivalence relation contained in both EQa and EQb. Both arguments must
    have the same length."""
    if len(EQa) != len(EQb): raise IndexMismatchError()
    
    N = len(EQa)
    alistsa = EQa.alists()
    alistsb = EQb.alists()
    components_b = dict([[root,bfs_component(alistsb,root)] for root in range(N) if EQb[root] < 0])
    ## components_b is a dictionary with items root : root/EQb, where root
    ## takes on all root values in EQb and root/EQb is a list of the elements
    ## in root's tree.
    
    component_list_r = []
    is_handled = [False] * N

    ## Basic procedure: 1: pop an element from the main queue.
    ## 2: find its EQa component and its EQb component.
    ## 2a: the intersection of these two is a component of EQr.
    ## 3. exhaust the remaining elements of the EQa component,
    ## treating it as a queue.

    for x in range(N):
        if is_handled[x]:
            pass
        else:
            cmp_a = bfs_component(alistsa,x)

            while len(cmp_a) > 0:
                x1 = cmp_a[0]
                cmp_b = components_b[EQb.find(x1)]

                cmp_r = [y for y in cmp_a if y in cmp_b]
                ## cmp_r is the intersection of cmp_a and cmp_b
                component_list_r.append(cmp_r)
                for y in cmp_r:
                    is_handled[y] = True

                cmp_a = [y for y in cmp_a if y not in cmp_r]
                ## deletes all elements of cmp_r from cmp_a

    return eqrel(N,blocks = component_list_r)

Python's list comprehension scheme is right up my alley, as you've no doubt already seen. I do sometimes worry about construction like the assignment of cmp_r near the end; abstractly, that could be \( \mathcal{O}(n^2) \), but I think Python has ways of making comprehensions like that run faster than the naive algorithm. (If someone knows if that's true, I'd love to know details. Maybe use some kind of hash and only check clashing pairs?)

So that's it for the methods I planned to write when I sat down to program. But wandering around in the documentation, I found something really fun, which I couldn't resist including. But that's for next time.

Blogging learning Python: equivalence relations I

 (Continuation here)

So I've finally broken down and started climbing up the learning curve for Python, just like everyone has been telling me to do for a couple of years now. Last summer's debacles with Octave should have been the last straw, but it's amazing what one can't get oneself to do when one actually got the result (despite the algorithm's taking a week to finally fail to compute what I asked for) and can move over to writing a paper instead of laboriously constructing examples.