[SciPy-user] C extension to manipulate sparse lil matrix

Andy Fraser afraser@lanl....
Thu Jan 8 18:32:26 CST 2009


I want to move some time critical bits of code for hidden Markov
models from python to C.  I've written code that works and uses sparse
matrices.  Next, I want to implement the "backward" algorithm in C.
As an intermediate step, I've coded/prototyped the manipulations that
I want to do on the internals of the sparse matrices using python.  I'll
append that code at the end here.

Now, I am trying to figure out how to manipulate lil sparse matrices.
In particular calling such a matrix "SM", and supposing that "t" is
the index for a row, I want to assign new arrays to "SM.rows[t]" and
"SM.data[t]".

I would be grateful if someone posted C code that interchanged two
rows of a lil sparse matrix.  I think I could glean what I need from
that example.  Since I'm new to C extensions, I'd like to see type
checking and reference counting done right too.

The basic recursion for the backward algorithm is

beta[t-1] = beta[t] {op1} Py[t] {op2} gamma[t] {op3} ScS

where beta[t-1], beta[t], and Py[t] are vectors, gamma[t] is a scalar,
and ScS is a matrix, and {op1} is element-wise multiplication of two
vectors, {op2} is division of a vector by a scalar, and {op3} is a
vector matrix product.

Here is my python code for the backward algorithm with sparse matrices:

======================================================================

def backsteps(N, T, gamma, Py_data, Py_rows, ScS_data, ScS_indices, ScS_indptr,
              beta_data,beta_rows):
    """ To imitate and check C.backsteps for debugging."""
    last_rows = numpy.array(range(N),numpy.int32)
    last_data = numpy.ones(N,numpy.float64)
    for t in xrange(T-1,-1,-1):
        beta_data[t] = last_data
        beta_rows[t] = last_rows
        gamma_t = gamma[t]
        Pyt_rows = Py_rows[t]
        Pyt_data = Py_data[t]
        mul_rows = []
        mul_data = []
        j0 = 0
        for i in xrange(len(Pyt_rows)):
            I = Pyt_rows[i]
            for j in xrange(j0,len(last_rows)):
                J = last_rows[j]
                if J == I:
                    mul_rows.append(I)
                    mul_data.append(Pyt_data[i]*last_data[j]/gamma_t)
                    j0 = j+1
                    break
                if J>I:
                    if j>j0:
                        j0 = j-1
                    break
        prod = numpy.zeros(N)
        for i in xrange(len(mul_rows)):
            I = mul_rows[i]
            for j in xrange(ScS_indptr[I],ScS_indptr[I+1]):
                J = ScS_indices[j]
                prod[J] += ScS_data[j]*mul_data[i]
        M = 0
        for i in xrange(N):
            if prod[i] != 0:
                M += 1
        last_rows = numpy.empty(M,numpy.int32)
        last_data = numpy.empty(M,numpy.float64)
        j = 0
        for i in xrange(N):
            if prod[i] != 0:
                last_rows[j] = i
                last_data[j] = prod[i]
                j += 1
    return


More information about the SciPy-user mailing list