[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